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1.0  INTRODUCTION 

Respiratory  failure  from  acute  lung  injury,  now  termed  the  acute  respiratory  distress  syndrome 
(ARDS),  accounts  for  4  million  ICU  days  annually  in  the  U.S.,  and  is  associated  with  high 
mortality  (up  to  40%)  in  both  military  and  civilian  populations.  Survivors  may  also  have 
substantial  morbidity,  with  long-term  physical  and  mental  health  impairments.  ARDS  thus  places 
significant  burdens  on  military  and  public  health  resources.  Since  ventilation  distribution  in  ARDS 
is  governed  by  a  heterogeneous  distribution  of  regional  mechanics,  the  most  appropriate  distending 
pressure,  ventilation  frequency,  or  tidal  volume  for  one  lung  region  may  not  necessarily  be  the 
same  for  another,  even  in  the  same  patient.  This  may  result  in  large  portions  of  an  injured  lung 
being  simultaneously  under-ventilated  or  over-ventilated,  with  poor  ventilation-to-perfusion 
matching  and  increased  dead  space.  The  goal  of  this  research  project  is  to  develop  a  novel 
method  for  mechanical  ventilation,  ‘Multi-Frequency  Oscillatory  Ventilation’  (MFOV)  that 
optimizes  gas  exchange  in  ARDS,  while  preserving  lung  protective  ventilation.  We 
hypothesize  that  lung  function  and  gas  exchange  will  be  significantly  improved  if  small  volume 
oscillations  are  applied  at  multiple  frequencies  simultaneously,  rather  than  at  a  single  frequency, 
due  to  more  even  distribution  of  ventilation  to  different  lung  regions  in  accordance  with  local 
mechanical  properties.  The  purpose  of  this  research  project  is  to  design  candidate  MFOV 
waveforms  for  the  acutely  injured  lung,  using  structurally  explicit  computational  models  of  the 
mammalian  respiratory  system.  We  then  use  dynamic  computed  tomographic  (CT)  imaging  to 
determine  whether  MFOV  improves  ventilation  distribution  and  gas  exchange  in  a  porcine  model 
of  ARDS.  Regional  ventilation  is  quantified  using  Xenon- enhanced  CT  in  these  pigs,  while 
regional  mechanical  properties  are  assessed  using  image  registration  and  measures  of  lung 
impedance.  We  expect  our  results  to  be  ultimately  translatable  and  testable  in  human  clinical  trials, 
with  potential  to  reduce  morbidity  and  mortality  associated  with  ARDS  and  other  heterogeneous 
lung  diseases. 

2.0  KEYWORDS 

Acute  lung  injury 

Acute  respiratory  distress  syndrome 
Blast  lung  injury 
Combat-related  lung  injury 
High  frequency  oscillatory  ventilation 
Multi-frequency  oscillatory  ventilation 
Oleic  acid  lung  injury 

3.0  ACCOMPLISHMENTS 
3.1  Goals  of  Project 

The  overall  goal  of  this  research  project  is  to  develop  a  novel  method  for  mechanical  ventilation, 
“Multi-Frequency  Oscillatory  Ventilation”  (MFOV)  that  optimizes  gas  exchange  in  ARDS  and 
other  forms  of  combat-related  lung  injury,  while  preserving  lung  protective  ventilation.  Our 
hypothesize  that  lung  function  and  gas  exchange  will  be  significantly  improved  if  small  volume 
oscillations  are  applied  at  multiple  frequencies  simultaneously,  rather  than  at  a  single  frequency, 
due  to  more  even  distribution  of  ventilation  to  different  lung  regions  in  accordance  with  local 
mechanical  properties.  The  goals  of  our  Specific  Aims  are:  1)  to  design  candidate  MFOV 
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waveforms  for  the  acutely  injured  lung,  using  structurally  explicit  computational  models  of  the 
mammalian  respiratory  system;  and  2)  to  use  dynamic  CT  imaging  to  establish  that  regional 
heterogeneity  is  the  mechanism  by  which  MFOV  improves  ventilation  distribution  and  gas 
exchange  in  a  porcine  model  of  ARDS. 

3.2  Goal  Accomplishments 

During  the  first  year  of  this  award  we  have  made  substantial  progress  along  our  research  plan,  both 
in  the  development  of  computational  models  for  the  design  MFOV  waveforms  in  acutely  lung 
injury,  as  well  as  the  confirmatory  animal  experiments  to  assess  regional  ventilation  distribution 
and  strain  during  MFOV.  Our  progress  is  detailed  for  both  Specific  Aims  1  and  2  below. 

3.2.1  Specific  Aim  1 

With  regard  to  the  computational  modeling  of  Specific 
Aim  1,  we  completed  the  development  of  a  three- 
dimensional  computational  model  of  a  canine  lung  to 
simulate  distributed  gas  flow  and  CO2  elimination 
during  oscillatory  ventilation.  Consistent  with  our 
hypothesis,  our  simulations  from  this  model 
demonstrated  that  ventilation  distribution  in  a 
heterogeneous  lung  is  spatially  clustered  and  dependent 
on  oscillatory  frequency  (Figure  1).  Such  regional 
differences  in  gas  exchange,  as  a  function  of  frequency, 
further  confirms  our  notion  that  use  of  MFOV  is  ideally 
suited  for  the  heterogeneously  injured  lung.  These 
results  were  recently  published  in  the  Journal  of  Applied 
Physiology  (11). 

To  determine  whether  such  model  could  be  used  to 
design  candidate  MFOV  waveforms  for  the  acutely 
injured  lung,  we  developed  a  Monte  Carlo  optimization 
algorithm  for  tuning  the  spectral  content  of  broadband 
oscillatory  flow  waveforms  according  to  regional 
mechanical  properties,  in  order  to  minimize 
parenchymal  strain  heterogeneity  (Figure  2).  Optimal 
combinations  of  frequencies,  amplitudes,  and  phases  in 
our  MFOV  waveforms  were  determined  according  to  frequency-dependent  distributions  of 
ventilation  throughout  the  heterogeneous  periphery  of  our  model.  We  found  that  the  superposition 
of  multiple  simultaneous  frequencies  provided  more  uniform  ventilation  distribution  in  our  model 
compared  to  single  frequency  oscillatory  ventilation,  with  less  potential  for  ventilator- induced  lung 
injury  (VILI).  These  results  were  presented  at  the  2016  Annual  Meeting  of  the  Biomedical 
Engineering  Society  in  Minneapolis,  Minnesota  (6). 

For  consistency  with  our  animal  experiments  of  Specific  Aim  2  (below),  we  have  now  extended 
these  computational  modeling  and  optimization  techniques  to  the  unique  anatomy  of  the  porcine 
lung,  to  investigate  the  potential  for  minimizing  distributed  acinar  strains  and  strain  rates  during 


Healthy  Injured 


1  2  5  10  20  50  100 

Frequency  [Hz] 


Figure  1:  Optimal  oscillation  frequencies 
for  CO2  elimination  in  a  healthy  and 
injured  canine  lung  model.  The  regional 
distribution  of  optimal  frequencies  in  the 
models  are  clustered  and  heterogeneous. 
These  are  simulated  assuming  supine 
position  of  the  animal,  with  the 
gravitational  vector  pointing  out  of  the 
page,  toward  the  reader.  Modified  from 
Herrmann  et  al.  (11). 
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Figure  2:  Convergence  of  the  MFOV  minimization  algorithm  from  randomized  initial  guesses  (indicated 
by  grey  lines)  to  local  minima  (colored  spheres).  In  this  case,  the  optimal  MFOV  solution  (red  sphere)  for 
3  frequencies  (/i,  /2,  /s  shown  in  Hz)  yielded  27%  reduction  in  parenchymal  strain  heterogeneity 
compared  to  the  optimal  solution  for  only  one  frequency.  The  resulting  optimized  flow  waveform  over  one 
second  is  also  shown  (inset). 


oscillatory  ventilation  (Figure  3).  Again,  our  simulations  demonstrated  that  MFOV  waveforms  are 
superior  to  traditional  single-frequency  HFOV  for  minimizing  strain  and  strain  rate  in  both  healthy 
and  injured  pig  lungs.  More  importantly,  we  found  that  an  MFOV  waveform  consisting  of  uniform 
flow  amplitudes  across  frequency  (i.e.,  hyperbolically  decreasing  volume  amplitudes  with 
increasing  frequency)  achieved  an  appropriate  balance  between  minimizing  parenchymal  strain 
vs.  strain  rate  in  the  porcine  lung.  These  results  will  be  presented  in  poster  format  at  the  2017 
American  Society  of  Anesthesiologists  Meeting  in  Boston,  Massachusetts  (14). 


Finally,  we  applied  our  computational  techniques  to  simulate  the  effects  of  a  very  specific  form  of 
lung  injury  with  military  relevance:  primary  blast-lung  injury.  A  computational  model  of  a  human 
airway  network  was  generated  using  central  airways  obtained  from  segmented  X-ray  computed 
tomographic  scans,  and  algorithmically  generated  peripheral  airways  (M.H.  Tawhai  et  al,  J  Appl 
Physiol  97(6):2310-2321,  2004).  The  model  consisted  of  60,494  cylindrical  airway  segments,  with 
30,243  terminal  bronchi  subtended  by  viscoelastic  acini.  Distributed  mechanical  properties  of 
tissues  and  airways  were  simulated  to  represent  typical  healthy  and  blast-injured  lungs,  the  latter 
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Figure  3:  Computational  model  of  porcine  lung  (left  panel).  Optimal  MFOV  volnme  amplitndes  at  each 
freqnency  are  shown  for  varying  cost  function  emphasis  on  strain  vs.  strain  rate,  in  simulated  healthy 
conditions  (center  panel)  and  heterogeneous  lung  injnry  (right  panel). 


characterized  by  bilateral  derecruitment  and  increased  tissue  stiffness,  focusing  on  the  perihilar 
regions  of  the  model  (Figure  4).  We  found  that  delivered  gas 
flow  was  heterogeneously  distributed  in  the  blast-injured  lung 
during  both  conventional  mechanical  ventilation  and  high- 
frequency  oscillatory  ventilation.  During  conventional 
mechanical  ventilation,  flow  was  distributed  primarily 
according  to  local  tissue  stiffness.  At  higher  frequencies,  the 
distribution  of  flow  became  increasingly  heterogeneous  and 
frequency-dependent,  with  some  regions  being  under- ventilated 
while  other  regions  experienced  substantially  greater  distension. 

These  results  were  presented  in  poster  form  at  the  2017  Military 
Health  System  Research  Symposium  in  Kissimmee,  Florida  (8). 

Ultimately,  our  unique  modeling  and  optimization  approaches 
will  allow  for  the  selection  of  patient- specific  MFOV 
waveforms,  especially  when  combined  with  experimental 
evidence  (Specific  Aim  2)  to  justify  physiologically-relevant 
emphasis  on  strain  vs.  strain  rate  to  minimize  risk  for  VILI. 
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Figure  4:  Healthy  human  airway 
tree,  along  with  three  levels  of 
blast  injury,  modeled  as 
increasing  levels  of  parenchymal 
tissue  stiffness  and  acinar 
derecruitment. 


Page  7 


3.2.2  Specific  Aim  2 

With  regard  to  Specific  Aim  2,  we 
have  completed  seven  animal 
experiments  to  date,  comparing 
conventional  mechanical  ventilation 
(CMV),  traditional  single-frequency 
HFOV,  and  MFOV  in  anesthetized 
and  paralyzed  pigs  before  and  after 
acute  lung  injury  across  several 
physiologic,  imaging,  and  gas 
exchange  metrics.  All  experiments 
have  been  approved  by  the  University 
of  Iowa  Institutional  Animal  Care  and 
Use  Committee  (Protocol  number 
5031314).  Measurements  of 
oscillatory  input  impedance  (13)  are 
obtained  under  baseline  conditions 
and  immediately  following  lung 
injury.  All  pigs  receive:  1)  CMV  at  a 
rate  of  20  min'^  and  tidal  volume 
between  10-12  mL  kg'^;  2)  HFOV 
delivered  at  5  Hz;  and  3)  MFOV 
delivered  using  uniform  flow 
amplitudes  at  5,  10,  15,  and  20  Hz. 

Each  ventilatory  /  oscillatory  modality 
is  applied  in  random  order,  under 
baseline  conditions  and  following 
lung  injury.  At  the  end  of  each  modality  period,  arterial  blood  gases  as  drawn,  and  CT  scans  are 
obtained  using  a  Siemens  Somatom  Force  scanner.  Figure  5  shows  example  transverse,  coronal, 
and  sagittal  CT  images  in  a  representative  pig  before  and  after  oleic  acid-induced  lung  injury. 
Similar  to  clinical  ARDS,  this  injury  model  is  characterized  by  patchy  alveolar  edema,  along  with 
generalized  increased  CT  density  that  is  greatest  in  the  dependent  lung  regions. 

During  the  course  of  these  experiments,  our  group  developed  a  revolutionary  dynamic  CT  image 
reconstruction  algorithm  to  assess  intratidal  derecruitment  and  overdistention  during  CMV  and 
MFOV  in  these  animals  with  acute  lung  injury.  This  frequency- selective  CT  (FSCT) 
reconstruction  can  identify  and  track  regions  of  the  lung  at  high  risk  for  additional  injury  due  to 
atelectrauma  or  volutrauma,  with  unparalleled  spatiotemporal  resolution.  Figure  6  shows  variation 
in  regional  aeration  over  the  course  of  a  positive  pressure  CMV  breath.  Our  technique  has  just 
recently  been  published  in  IEEE  Transactions  in  Medical  Imaging  (4).  In  collaboration  with  Drs. 
Gary  E.  Christensen  and  Joseph  M.  Reinhardt  at  the  University  of  Iowa,  we  have  applied  a  novel 
4-D  image  registration  technique  to  our  unique  ESCT  image  sequences,  to  reconstruct  temporally- 
resolved  structural  deformations  during  periodic  ventilation,  along  with  regional  maps  of 
parenchymal  expansion  and  contraction  (15).  Periodic  expansion  and  contraction  are  assessed  by 
the  discrete  Eourier  transform  of  the  normalized  time- varying  Jacobian  determinants.  Eigure  7 
shows  the  distribution  of  strain  amplitudes  in  a  representative  pig  during  HEOV  delivered  at  5  Hz, 


Baseline  Oleic  Acid  Injury 


Figure  5:  CT  images  in  a  representative  pigs  before  (left  panels) 
and  after  (right  panels)  lung  injury,  shown  in  transverse  (top 
row),  coronal  (middle  row),  and  sagittal  (bottom  row)  sections. 
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Figure  6:  Changes  in  regional  aeration  during  porcine  oleic  acid  lung  injury  using  dynamic  FSCT 
reconstruction.  Ranges  of  CT  Hounsfield  density  correspond  to  levels  of  lung  aeration.  Transverse  sections 
of  the  FSCT  reconstruction  are  shown  at  four  time  points  during  a  positive  pressure  CMV  hreath.  Note 
that  26  %  of  the  lung  cyclically  alters  between  normally  and  poorly  aerated. 


as  well  as  during  MFOV.  These  distributions  tend  to  vary  primarily  in  the  ventral-dorsal  direction, 
in  accordance  with  the  gravitational  field  in  supine  position.  Regional  strain  amplitudes  during 
MFOV  varies  with  frequency,  with  different  lung  regions  selectively  filtering  the  harmonic 
frequency  content  of  the  broadband  oscillatory  flow.  Of  note,  the  total  volumetric  strain  (peak-to- 
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Figure  7:  Regional  distrihutions  of  volumetric  strain  amplitudes  during  HFOV  and  MFOV  shown  in 
transverse  cross-section,  using  discrete  Fourier  transform  of  normahzed  time-varying  Jacohian 
determinants. 


Page  9 


PIG-06  PIG-07  PIG-08  PIG-09  PIG-10  PIG-11 


Figure  8:  Distributions  of  xenon  wash-in  time  constants  in  six  pigs  under  baseline  conditions  during 
mechanical  ventilation  using  either  CMV,  HFOV,  or  MFOV.  Nonlinear  regression  was  used  to  fit 
parameters  of  exponential  saturation  to  respiratory-gated  CT  image  sequences 

peak  mean  ±  standard  deviation)  during  HFOV  (7.9%  ±  3.1%)  and  MFOV  (6.9%  ±  2.5%)  were 
substantially  lower  than  during  CMV  (31.3%  +/-  11.6%).  Our  preliminary  results  indicate  that 
parenchymal  stretch  during  oscillatory  ventilation  is  regionally  heterogeneous  and  frequency- 
dependent.  Moreover,  the  broadband  spectral  content  associated  with  MFOV  may  enhance  gas 
transport  in  the  presence  of  periodic  parenchymal  deformation  with  higher  harmonics.  Such 
behavior  may  be  further  adjusted  to  compensate  for  patient- specific  regional  heterogeneity.  These 
results  will  be  presented  at  the  2017  meeting  of  the  Biomedical  Engineering  Society  in  Phoenix, 
Arizona  (7). 

In  conjunction  with  the  FSCT  imaging  described  above,  we  also  acquired  timed  sequences  of  CT 
images  during  xenon  wash-in  and  wash-out  in  each  pig  to  assess  regional  ventilation  distribution. 
Image  acquisition  was  gated  to  the  periodic  mechanical  ventilation  cycle  during  CMV,  HFOV,  or 
MFOV.  All  images  in  each  sequence  were  pairwise  registered  to  their  respective  initial  images,  to 
compensate  for  non-respiratory  sources  of  motion  such  as  cardiogenic  oscillations.  Regions  of 
interest  were  defined  as  cubes  spanning  four  voxels  on  each  side.  Within  each  region  of  interest, 
nonlinear  regression  was  used  to  fit  exponential  saturation  or  decay  curves  to  a  time  series  of 
spatially- averaged  voxel  intensities  during  xenon  wash-in  or  wash-out,  respectively.  Figure  8 
shows  representative  slices  from  xenon  wash-in  image  analysis  performed  in  six  pigs  under 
baseline  conditions.  These  results  demonstrate  the  feasibility  of  our  XeCT  imaging  protocol  in 
healthy  subjects,  and  furthermore  suggest  noticeable  qualitative  differences  in  the  distribution  of 
xenon  equilibration  time  constants  between  the  ventilation  modalities.  In  particular,  gas  transport 
during  MFOV  appears  to  be  least  affected  by  gravitational  forces  compared  to  CMV  or  HFOV. 
Additionally,  HFOV  tends  to  produce  larger  time  constants  than  CMV  or  MFOV,  suggesting 
overall  slower  rates  of  xenon  equilibration  during  HFOV. 


Page  10 


3.3  Opportunities  for  Training  and  Professional  Development 

Jacob  Herrmann,  M.S.,  is  currently  a  Ph.D.  candidate  in  the  Department  of  Biomedical 
Engineering  at  the  University  of  Iowa,  and  has  been  working  under  the  direct  supervision  of  Dr. 
Kaczka  on  all  aspects  of  this  project.  Mr.  Herrmann  has  assisted  with  all  animal  experiments  to 
date,  and  has  been  managing  the  analyses  for  the  CT  imaging  and  impedance  data  in  pigs.  As  a 
part  of  his  Ph.D.  research,  he  has  developed  the  dynamic  FSCT  image  reconstruction  algorithm  to 
quantify  lung  parenchymal  motion  during  mechanical  ventilation  and  oscillation  (4),  as  well  as  4- 
D  image  registration  to  assess  mechanical  strain  heterogeneity  in  the  injured  porcine  lungs.  Along 
with  our  consultant  Dr.  Merryn  Tawhai  of  the  University  of  Auckland,  Mr.  Herrmann  has 
developed  a  three-dimensional  computational  model  of  a  porcine  lung  with  high  anatomic  fidelity, 
based  on  a  central  airway  tree  segmented  from  a  CT  image  of  a  supine  pig.  Simulations  of 
ventilation  distribution  using  this  model  has  been  crucial  to  our  understanding  how  MFOV 
enhances  gas  exchange  in  the  injured  lung,  and  how  spectral  content  of  MFOV  waveforms  can  be 
further  optimized  according  to  the  relative  contributions  of  acinar  strain  vs.  strain  rate.  During  the 
first  year  of  this  grant,  Mr.  Herrmann  has  been  the  first  author  of  two  publications  (4,  11),  and  has 
another  manuscript  in  press  (10). 

3.4  How  were  the  results  disseminated  to  the  communities  of  interest? 

We  have  published  one  article  (11)  and  have  another  article  (10)  accepted  for  publication  in  the 
Journal  of  Applied  Physiology.  The  audience  for  this  journal  consists  mostly  of  organ-level 
physiologists  and  biomedical  engineers.  Our  paper  focusing  on  the  FSCT  technique  (4)  was 
published  in  IEEE  Transactions  in  Medical  Imaging,  with  an  audience  of  clinical  radiologists, 
imaging  scientists,  and  medical  physicists.  Our  work  from  this  contract  has  also  been  presented  at 
meetings  of  the  Biomedical  Engineering  Society  (1,7,  9),  the  American  Thoracic  Society  (3),  the 
Association  of  University  Anesthesiologists  (5),  the  International  Flow-Volume  Underworld 
Meeting  (2,  12)  and  the  Military  Health  System  Research  Symposium  (8). 

3.5  What  do  you  plan  to  do  during  the  next  reporting  period  to  accomplish  the  goals? 

We  plan  to  finish  the  experimental  protocols  in  the  three  remaining  pigs.  We  are  now  in  the  in  the 
ongoing  process  of  analyzing  and  comparing  all  gas  exchange  data,  physiological  variables  and 
indices,  as  well  as  dynamic  CT  imaging  data.  A  very  time  consuming  process  is  the  segmentation 
of  CT  images  of  the  injured  pig  lungs,  since  the  Hounsfield  density  of  parenchymal  consolidation 
and  edema  fluid  are  difficult  to  distinguish  from  other  chest  wall  structures.  This  makes  semi¬ 
automatic  segmentation  algorithms  ineffective,  and  we  must  rely  on  manual  segmentation. 
However,  we  expect  to  have  this  completed  by  the  closeout  of  the  final  progress  report. 

4.0  IMPACT 

4.1  What  was  the  impact  on  the  development  of  the  principal  discipline  of  the  project? 

MFOV  a  promising  and  innovative  approach  to  lung  protective  ventilation,  and  its  use  in  patients 
will  require  considerable  work  to  understand  its  application  in  pathophysiologies  relevant  to 
ARDS  and  other  forms  combat-related  lung  injury.  The  mechanical  heterogeneity  of  the  injured 
lung  has  important  implications  for  optimal  ventilation  protocols  and  treatment  strategies.  The 
novel  modeling,  experimental,  and  imaging  techniques  detailed  in  this  progress  report  address 
fundamental  questions  regarding  the  mechanism  by  which  MFOV  improves  gas  exchange 
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compared  to  CMV  or  traditional  HFOV.  The  immediate,  short-term  impact  of  this  project  is  the 
development  of  MFOV  waveforms  that  have  far  more  potential  to  reduce  the  risk  of  ventilator- 
associated  lung  injury  compared  to  CMV,  as  well  as  a  mechanistic  understanding  of  gas  transport 
during  oscillatory  ventilation.  Our  work  also  provides  a  requisite,  scientific  basis  for  the  eventual 
use  of  MFOV  in  patients,  since  our  results  are  easily  translatable  and  testable  in  human  clinical 
trials.  MFOV  will  have  long-term  impact  on  guiding  therapeutic  and  technological  developments 
for  future  research  in  ARDS,  with  potential  to  significantly  impact  the  morbidity  and  mortality 
associated  with  this  deadly  syndrome. 

4.2  What  was  the  impact  on  other  disciplines? 

As  detailed  in  Section  3.0  of  this  progress  report,  our  new  frequency-selective  CT  (FSCT)  dynamic 
image  reconstruction  algorithm  enables  imaging  of  both  low-  and  high-frequency  dynamic 
periodic  motions  at  specified  frequencies,  with  minimal  extraneous  motion  artifact.  This  FSCT 
approach  offers  improved  imaging  fidelity  of  dynamic  subject  motion  compared  to  more 
conventional  prospective  or  retrospective  gating  approaches.  The  characterization  of  dynamic 
deformations  of  thoracic  structures  will  have  enormous  impact  in  the  fields  of  clinical  radiology 
and  cardiac  imaging.  When  combined  with  image  registration  techniques,  FSCT  may  provide 
detailed  four-dimensional  information  of  distributed  tissue  deformation  of  anatomic  structures 
previously  inaccessible  using  standard  reconstruction  approaches. 

4.3  What  was  the  impact  on  technology  transfer? 

Nothing  as  yet  to  report 

4.4  What  was  the  impact  on  society  heyond  science  and  technology? 

ARDS  has  a  major  impact  in  the  United  States,  with  an  estimated  190,000  cases  and  74,000  deaths 
annually.  Significant  reductions  in  mortality  have  been  realized  by  the  use  of  lung  protective 
mechanical  ventilation  protocols  in  which  PEEP  is  used  to  recruit  the  lung  and  prevent  repetitive 
end-expiratory  opening  and  closing  of  airspaces,  and  smaller  tidal  volumes  prevent  end- inspiratory 
over-distension.  MFOV  represents  a  promising  and  innovative  approach  to  lung  protective 
ventilation,  although  its  use  in  patients  will  require  considerably  more  work  to  understand  its 
application  in  pathophysiologies  relevant  to  ARDS  and  combat-related  lung  injury.  This  project 
will  ultimately  provide  a  solid,  scientific  basis  for  the  rational  use  of  this  MFOV  in  critically  ill 
patients.  The  thrust  of  this  DoD  CDMRP  contract  is  to  demonstrate  that  MEOV  results  in  improved 
gas  exchange  while  minimizing  additional  injury  in  an  animal  model  of  combat-related  lung  injury. 
However,  MFOV  may  have  more  far-reaching  implications  for  both  pulmonary  medicine  and 
anesthesia.  Eor  example,  MEOV  may  not  be  limited  to  a  treatment  solely  for  ARDS,  but  may  be 
utilized  in  the  ventilator  management  of  other  heterogeneous  lung  diseases,  such  as  asthma, 
COPD,  or  pneumonia.  The  possibility  that  MEOV  can  more  efficiently  penetrate  ‘difficult-to- 
reach’  regions  of  the  lung  also  has  implications  for  the  optimal  delivery  of  aerosols  and  drugs, 
such  as  beta  agonists,  steroids,  or  even  inhaled  volatile  anesthetics. 
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5.0  CHANGES  /  PROBLEMS 

5.1  Changes  in  approach  and  reasons  for  change 

Our  original  experimental  protocol  indicated  that  four  ventilatory  /  oscillatory  modalities  would 
be  used  in  each  animal:  1)  CMV;  2)  single-frequency  HFOV;  3)  a  generic  MFOV  consisting  of 
uniform  flow  amplitudes  (or  hyperbolically-decreasing  volume  amplitudes);  and  4)  a  custom- 
designed  MFOV  waveform  based  on  the  modeling  and  simulation  results  of  Specific  Aim  1. 
However,  our  computer  modeling  results  indicated  that  an  MFOV  waveform  consisting  of 
hyperbolically-decreasing  volume  amplitudes  provided  a  near-optimal  solution  for  most  levels  of 
lung  injury  (Figure  3).  Thus,  we  were  able  to  simplify  our  experimental  protocol  using  only  three 
modalities,  since  our  original  generic  MFOV  waveform  was  found  to  be  optimal.  This  was  a 
relatively  minor  change  to  our  initial  experimental  design. 

5.2  Actual  or  anticipated  problems  or  delays  and  actions  or  plans  to  resolve  them 

Nothing  to  report 

5.3  Changes  that  had  a  significant  impact  on  expenditures 

Nothing  -  we  are  within  our  budget,  despite  doing  more  experiments  than  planned  in  the  first  year. 
As  detailed  in  Section  5.1,  the  experimental  protocol  was  simplified  from  comparing  four 
ventilatory  modalities  to  just  three,  since  our  ‘generic’  MFOV  waveform  yielded  optimal  solutions 
based  on  our  computer  modeling.  Thus  the  duration  for  each  experiment  was  shortened,  and  thus 
reduced  expenditures.  We  budgeted  for  3  experiments  in  Year  1  and  7  experiments  in  Year  2, 
however,  we  were  able  to  get  completed  protocols  in  7  experiments  in  Year  1  with  sufficient 
budget  for  4-5  more  experiments  in  Year  2.  Analysis  of  the  completed  experiments  from  Year  1 
is  continuing  into  Year  2. 


Budget  category 

Total  budget 

Cumulative 
Expenditures  for 
Year  1 

Remaining  for 

Year  2 

Animal  experiment 
costs  (7  experiments) 

$76,982.50 

$55,642.45 

$21,340.05 

5.4  Significant  changes  in  the  use  or  care  of  human  subjects,  vertebrate  animals,  biohazards, 
and/or  select  agents 

Nothing  to  report. 

6.0  PRODUCTS 

6.1  Journal  Publications 

Herrmann  J,  Tawhai  MH,  and  Kaczka  DW.  Regional  gas  transport  in  the  heterogeneous  lung 
during  oscillatory  ventilation.  Journal  of  Applied  Physiology  121:  1306-1318,  2016; 
Acknowledgement  of  federal  support:  yes 
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Herrmann  J,  Hoffman  EA,  and  Kaczka  DW.  Frequency-Selective  Computed  Tomography: 
Applications  During  Periodic  Thoracic  Motion.  IEEE  Trans  Med  Imaging  36:  1722-1732,  2017; 
Acknowledgement  of  federal  support:  yes 

Herrmann  J,  Tawhai  MH,  and  Kaczka  DW.  Parenchymal  strain  heterogeneity  during  oscillatory 
ventilation:  Why  two  frequencies  are  better  than  one.  Journal  of  Applied  Physiology  (in  press); 
Acknowledgement  of  federal  support:  yes 

6.2  Books  or  other  non-periodical,  one-time  publications 

Nothing  to  report 

6.3  Other  Publications,  Conference  Papers,  and  Presentations 

Herrmann  J,  and  Kaczka  DW.  Optimization  of  Spectral  Content  in  Oscillatory  Ventilator 
Waveforms.  In:  2016  Annual  Meeting  of  the  Biomedical  Engineering  Society.  Minneapolis,  MN: 
2016.  Acknowledgement  of  federal  support:  yes 

Herrmann  J.  What  can  frequency-selective  CT  imaging  show  us  about  regional  lung  deformation? 
Elow -Volume  Underworld,  Boston,  MA,  2017.  Acknowledgement  of  federal  support:  yes 

Kaczka  DW.  Parenchymal  Strain  During  Oscillatory  Ventilation:  Why  Two  Frequencies  are  Better 
than  One.  Elow-Volume  Underworld,  Boston,  MA,  2017.  Acknowledgement  of  federal  support: 
yes 

Herrmann  J,  and  Kaczka  DW.  Minimizing  Parenchymal  Strain  Heterogeneity  During  Oscillatory 
Ventilation.  In:  Association  of  University  Anesthesiologists  64th  Annual  Meeting.  Washington, 
D.C.:  2017.  Acknowledgement  of  federal  support:  yes 

Herrmann  J,  Fuld  MK,  Hoffman  EA,  and  Kaczka  DW.  Temporally  Resolved  Computed 
Tomographic  (CT)  Imaging  to  Assess  Regional  Parenchymal  Strain  during  Oscillatory 
Ventilation.  2017  Meeting  of  the  American  Thoracic  Society.  American  Journal  of  Respiratory 
and  Critical  Care  Medicine  2017;  195:A6522.  Acknowledgement  of  federal  support:  yes 

Herrmann  J,  Shao  W,  Reinhardt  J,  Hoffman  E,  Christensen  G,  and  Kaczka  DW.  Frequency- 
Selective  CT  Image  Registration  for  Assessment  of  Regional  Periodic  Fung  Deformation.  2017 
Meeting  of  the  Biomedical  Engineering  Society,  Phoenix  AZ,  Acknowledgement  of  federal 
support:  yes 

Herrmann  J,  Tawhai  M,  and  Kaczka  DW.  Multi- Objective  Optimization  of  Multi-Frequency 
Oscillatory  Ventilation.  2017  Meeting  of  the  Biomedical  Engineering  Society,  Phoenix  AZ, 
Acknowledgement  of  federal  support:  yes 

Hadjarevic  B,  Herrmann  J,  Kaczka  DW.  Modeling  the  Frequency  Response  of  a  Proportional 
Solenoid  Valve  for  Oscillatory  Ventilation.  2017  Meeting  of  the  Biomedical  Engineering  Society, 
Phoenix  AZ,  Acknowledgement  of  federal  support:  yes 
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6.4  Website  or  other  Internet  site 


https://www.researchgate.net/project/Optimizing-Ventilation-Distribution-and-Gas-Exchange-in- 
Combat- Related- Lung- Injury-Using-Multi-Frequency-Oscillation 

6.5  Technologies  or  techniques 

Nothing  to  report 

6.6  Inventions,  patent  applications  and/or  licenses 

Nothing  to  report 

6.7  Other  Products 

Nothing  to  Report 

7.0  PARTICIPANTS  &  OTHER  COLLABORATING  ORGANIZATIONS 
7.1  What  individuals  have  worked  on  the  project? 


Name: 

David  W.  Kaczka,  M.D.,  Ph.D. 

Project  Role: 

Principal  Investigator 

Researcher  Identifier  (ORCID  ID): 

0000-0003-4378-5242 

Nearest  person  month  worked: 

3.2 

Contribution  to  Project: 

Dr.  Kaczka  is  responsible  for  the  overall 
direction  of  the  project  as  well  as  the 
experimental  measurements  of  lung 
impedance  and  CT  imaging.  Dr.  Kaczka  is 
responsible  for  oversight  of  all  publications 
and  presentations  related  to  the  project.  Dr. 
Kaczka  supervised  and  participated  in  the 
animal  experiments  with  Mr.  Jacob  Herrmann 
(Ph.D.  student). 

Funding  Support: 

Nm 

Name: 

Jacob  Herrmann,  M.S. 

Project  Role: 

Graduate  Student 

Researcher  Identifier  (ORCID  ID): 

0000-0001-5046-5592 

Nearest  person  month  worked: 

12 

Contribution  to  Project: 

Mr.  Herrmann  has  performed  the 
computational  modeling  for  assessment  of  gas 
transport  and  exchange  in  porcine  lungs  in 
collaboration  with  Dr.  Kaczka  and  Dr.  Tawhai. 
He  has  also  assisted  in  the  animal  experiments 
and  analysis  of  the  imaging  and  impedance 
data. 

Funding  Support: 

University  of  Iowa,  Department  of  Anesthesia 
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Name: 

Merryn  Tawhai,  Ph.D. 

Project  Role: 

Consultant 

Researcher  Identifier  (ORCID  ID): 

0000-0002-3211-6337 

Nearest  person  month  worked: 

160  hours  over  Year  1  (1.0  person  month) 

Contribution  to  Project: 

Advisor  for  modeling  projects  in  Year  1, 
completed 

Funding  Support: 

NIH,  Royal  Society  of  New  Zealand,  Health 
Research  Council  of  New  Zealand,  The 
Tertiary  Education  Commission  of  New 
Zealand 

Name: 

Eric  Hoffman,  Ph.D. 

Project  Role: 

Co-Investigator 

Researcher  Identifier  (ORCID  ID): 

0000-0001-8456-9437 

Nearest  person  month  worked: 

0.6 

Contribution  to  Project: 

Dr.  Hoffman  closely  collaborated  in  all  aspects 
of  the  computational  modeling,  animal 
experiments,  CT  scanning  and  image 
processing 

Funding  Support: 

NIH 

7.2  Has  there  been  a  change  in  the  active  other  support  of  the  PI  or  key  personnel  since  the 
last  reporting  period? 

The  PI,  David  Kaczka,  MD,  PhD,  has  added  effort  on  two  projects,  however  overall  research  effort 
is  unchanged  from  60%  research  effort,  40%  clinical  effort.  There  is  no  overlap  between  the 
current  project  and  these  projects. 

2015-2019  NIH/NHLBI 1R01HL126838 

“Defining  the  Role  of  Early  Pulmonary  Vascular  Disease  in  COPD” 

Role:  Co-Investigator  10%  FTE 

This  project  evaluates  the  role  of  imaging-based  measures  of  pulmonary 
vascular  structure  and  function  to  identify  smokers  at  risk  of  rapid 
emphysema  progression  and/or  at  risk  of  having  frequent  acute 
exacerbations  of  their  COPD. 

2012-2017  NIH  5R01HE1 12986 

“Multi-Center  Structural  &  Functional  Quantitative  CT  Pulmonary 
Phenotyping” 

Role:  Co-Investigator  10%  FTE 

This  bioengineering  research  partnership  takes  advantage  of  emerging 
acquisition  technique  of  multi-spectral  computed  tomography  (dual  energy 
CT:  DECT),  careful  evaluation  of  dose  lowering  methods,  and  novel 
approaches  to  statistical  cluster  analysis  to  expand  the  biomarkers  used  in 
multi-center  studies  to  identify  sub-populations  of  lung  disease. 
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The  Co-I,  Eric  Hoffman,  PhD,  has  modified  his  effort  on  various  NIH  projects.  There  is  no  overlap 
between  the  current  project  and  any  of  his  other  projects. 

7.3  What  other  organizations  were  involved  as  partners? 

Organization  Name:  University  of  Auckland 

Location  of  Organization:  Auckland,  New  Zealand 

Partner’s  contribution  to  the  project:  Consultant,  Merryn  Tawhai,  Ph.D. 

Dr.  Tawhai  advised  Dr.  Kaczka  and  Mr.  Herrmann  on  aspects  of  the  modeling  that  was  completed 
in  Year  1.  Her  commitment  was  as  expected  and  the  full  budget  of  $11,840  had  been  paid.  No 
further  effort  is  required  by  Dr.  Tawhai  on  this  project. 

8.0  SPECIAL  REPORTING  REQUIREMENTS 

Nothing  to  report 
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Modeling  the  Frequency  Response  of  a  Proportional  Solenoid  Valve  for  Oscillatory  Ventilation 

Bakir  Hajdarevic*,  Jacob  Herrmann*,  and  David  Kaczka* 

'University  of  Iowa,  Iowa  City,  lA 


Introduction:  Proportional  solenoid  (PSOL)  valves  are  an  important  component  of  modern  mechanical  ventilators***.  However 
the  use  of  PSOL  valves  in  high  frequency  oscillators  is  limited  due  to  incomplete  characterization  of  their  dynamic  responses,  as 
well  as  the  necessity  to  incorporate  them  into  sophisticated  closed-loop  systems  for  precise  airway  pressure  control  and  patient 
safety*^*.  The  goal  of  this  study  was  to  characterize  the  dynamic  response  of  high-flow  PSOL  valve  using  a  linear  transfer 
function,  which  may  allow  for  robust  simulation  and  design  of  closed-loop  controllers  for  various  ventilatory  and  oscillatory 
modalities  in  patients. 

Materials  and  Methods:  Measurements  were  made  in  an  ASCO  Posiflow™  PSOL  valve  (Model  SD8202G027V)  connected  in 
series  to  a  Bellofram  Type  70  Air  Regulator  (960-131-010).  Maximum  orifice  diameter  of  the  valve  was  7  mm,  allowing  for  a 
steady  peak  flow  delivery  of  4.7  L  sec  *.  Pulse  width  modulated  (PWM)  current  to  the  PSOL  coil  was  provided  by  ASCO’s 
Electronic  Control  Unit  (ECU  Model  8908A001)  at  a  frequency  of  300  Hz.  A  digital-to-analog  (D/A)  converter  (USB-6008, 
National  Instruments)  was  used  to  generate  bandlimited  white  noise  (0.01-40  Hz),  as  well  as  three  different  pseudorandom  signals 
consisting  of  discrete  sinusoids  at  mutually  prime  frequencies  (0.195-37.4  Hz).  The  output  of  the  D/A  converter  was  low  pass 
filtered  at  40  Hz  (8-pole  Butterworth),  and  presented  as  the  input  control  voltage  to  the  ECU.  Output  flow  from  the  PSOL  valve 
was  measured  with  a  pneumotachograph  (Hans  Rudolph  4700A)  coupled  to  a  0-2  cniH20  differential  pressure  transducer 
(Celesco  LCVR-0002).  Both  the  input  voltage  and  output  flow  waveforms  were  low  pass  filtered  at  40  Hz,  and  sampled  by  an 
analog-to-digital  converted  at  100  Hz.  The  voltage-flow  transfer  function  of  the  combined  ECU-PSOL  system  was  determined 
using  an  overlap-average  periodogram  technique,  which  was  then  fitted  to  a  generalized  parametric  transfer  function  of  the  form: 

M  /  N 

r(^)=/^n(^+aj/n(^+i3j  (1) 

m=l  /  «=1 

where  a,„  and  p„  denote  the  zeroes  and  poles,  respectively,  of  T{s),  and  M  <  N.  The  parameters  K,  a,„,  and  p„  were  estimated 
using  a  nonlinear  gradient  search  technique  (MATLAB  v9.0.  The  Mathworks  Inc.).  The  optimal  number  of  zeros  and  poles  for 
the  ECU-PSOL  system  (M,  N)  were  determined  based  on  the  Akaike  information  criterion  (AIC). 

Results  and  Discussion:  The  PSOL  valve  exhibited  a  steady-flow 
hysteresis  of  -18%,  while  its  dynamic  magnitude  and  phase  responses  were 
flat  out  to  10  Hz  (Eigure  1).  Broadband  measurement  of  total  harmonic  flow 
distortion*^*  was  approximately  11%.  The  frequency  response  of  the  ECU- 
PSOL  system  was  adequately  described  by  Equation  1  with  M  =  3  and  N  = 

5. 


Conclusions:  The  dynamic  frequency  response  of  a  high-flow  PWM  PSOL 
valve  can  be  adequately  described  using  a  linear  polynomial  transfer 
function.  Such  a  characteristic  will  allow  for  robust  simulation  and  design 
of  a  closed-loop  controller  in  a  device  capable  of  generating  various 
ventilatory  and  oscillatory  modalities  in  patients  with  acute  respiratory 
failure. 
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Figure  1  Frequency  response  of 
PSOL  valve  driven  with  white  noise 
(gray  symbols)  and  pseudorandom 
signal  (hlack  symbols)  inputs,  along 
with  corresponding  fit  of  Equation  1 
to  pseudorandom  data  (solid  lines). 
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Abstract  Body 
Rationale: 

We  have  recently  proposed  the  use  of  multi-frequency  oscillatory  ventilation  (MFOV)  as  a  modification  of  high  frequency 
oscillatory  ventilation  (FI FOV)  for  improving  ventilation  distribution  and  minimizing  ventilator-induced  lung  injury  (Kaczka  et 
al.,  Anesf/ies/o/ogy  123(6):  1394-1 403,  2015).  Direct  observation  of  dynamic  regional  strain  during  MFOV  is  not  possible 
using  conventional  X-ray  CT  imaging,  due  to  motion  artifact  blurring  and  extraneous  cardiogenic  disturbances.  These 
artifacts  are  especially  problematic  during  oscillation  at  frequencies  greater  than  the  rotation  frequency  of  the  X-ray  source 
and  detector  array.  In  this  study,  we  developed  a  reconstruction  algorithm  for  spatial  and  temporal  resolution  of  periodic 
lung  motion  during  CT  scanning,  enabling  quantification  of  regional  parenchymal  strain  during  FIFOV  or  MFOV. 

Methods: 

Sedated  and  paralyzed  pigs  received  mechanical  ventilation  (conventional,  FIFOV,  or  MFOV)  using  ventilatory  frequencies 
between  0.4  to  20  Hz  (i.e.,  24  to  1 200  min‘^ ).  CT  projections  were  acquired  continuously  for  30  to  1 00  seconds  using  a 
SOMATOM  Force  scanner  (Siemens  Healthineers,  Forchhiem,  Germany)  rotating  at  4  Hz,  using  80  kVp  and  150  mAxs  (total 
dose  per  acquisition  between  345  to  1152  mGy).  CT  projections  were  reorganized  according  to  oscillatory  frequency  before 
cone-beam  reconstruction,  yielding  a  series  of  volumetric  images  (isotropic  0.3  mm  voxel  size,  192  slices,  240  mm  field-of- 
view)  at  10  distinct  time-points  during  periodic  motion,  with  minimal  extraneous  motion  artifact.  Piecewise  four-dimensional 
image  registration  was  used  to  determine  spatial  parenchymal  deformation  between  consecutive  time-points. 

Results: 

Our  reconstruction  algorithm  produced  image  sequences  of  clearly  defined  lung  structure  at  distinct  time  points  throughout 
the  oscillatory  ventilation  cycle  (Figure  1).  Substantial  reductions  in  reconstruction  error  were  achieved  for  the  longest 
scanning  durations.  Cardiogenic  oscillations  resulted  in  minor  motion  artifact  due  to  the  oscillatory  frequency  not  being  an 
integer  multiple  of  heart  rate.  Image  registration  revealed  heterogeneous  distributions  of  parenchymal  strain  magnitudes 
and  phases  during  oscillatory  ventilation. 

Conclusions: 

We  have  developed  a  novel  method  for  reconstructing  temporally  resolved  images  of  oscillatory  lung  deformation  during 
continuous  X-ray  CT  acquisition,  with  minimal  distortion  from  extraneous  motion.  Cur  method,  coupled  with  image 
registration  techniques,  enables  direct,  non-invasive  quantification  of  regional  parenchymal  strain  in  vivo  during  high 
frequency  and  multi-frequency  oscillatory  ventilation. 
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Applications  During  Periodic  Thoracic  Motion 

Jacob  Herrmann,  Student  Member,  IEEE,  Eric  A.  Hoffman,  Member,  IEEE, 
and  David  W.  Kaczka*,  Member,  IEEE 


Abstract — We  seek  to  use  computed  tomography  (CT)  to 
characterize  regional  lung  parenchymal  deformation  during 
high-frequency  and  multi-frequency  oscillatory  ventilation. 
Periodic  motion  of  thoracic  structures  results  in  artifacts 
of  CT  images  obtained  by  standard  reconstruction  algo¬ 
rithms,  especially  for  frequencies  exceeding  that  of  the 
X-ray  source  rotation.  In  this  paper,  we  propose  an  acquisi¬ 
tion  and  reconstruction  technique  for  high-resolution  imag¬ 
ing  of  the  thorax  during  periodic  motion.  Our  technique 
relies  on  phase-binning  projections  according  to  the  fre¬ 
quency  of  subject  motion  relative  to  the  scanner  rotation, 
prior  to  volumetric  reconstruction.  The  mathematical  theory 
and  limitations  of  the  proposed  technique  are  presented, 
and  then  validated  in  a  simulated  phantom  as  well  as  a 
living  porcine  subject  during  oscillatory  ventilation.  The 
4-D  image  sequences  obtained  using  this  frequency- 
selective  reconstruction  technique  yielded  high-spatio- 
temporal  resolution  of  the  thorax  during  periodic  motion. 
We  conclude  that  the  frequency-based  selection  of  CT  pro¬ 
jections  is  ideal  for  characterizing  dynamic  deformations  of 
thoracic  structures  that  are  ordinarily  obscured  by  motion 
artifact  using  conventional  reconstruction  techniques. 

Index  Terms — Heart,  image  acquisition,  image  recon¬ 
struction  -  analytical  methods,  lung,  motion  compensation 
and  analysis,  tracking  (time  series  analysis).  X-ray  imaging 
and  computed  tomography. 


Glossary  of  Nomenclature 


FSCT 

Frequency-selective  computed  tomography 

/rot 

Frequency  of  scanner  rotation. 

/sub 

Frequency  of  subject  motion. 

k 

Index. 

^bin 

Index  of  subject  phase  bin. 

^pro 

Index  of  projection  angle. 

-^bin 

Number  of  subject  phase  bins. 

^pro 

Number  of  projections  acquired  per 
scanner  rotation. 
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Nrot  Total  scanning  duration,  measured  in  number  of 
scanner  rotations. 

p  Index  of  sequentially  acquired  projections. 

Rpro  Number  of  projections  acquired  before  repetitious 
sampling. 

Rrot  Number  of  scanner  rotations  completed  before 
repetitious  sampling. 

Rsub  Number  of  subject  motion  cycle  completed  before 
repetitious  sampling. 
t  Time. 

tic  Times  corresponding  to  the  same  subject  phase. 
tp  Time  of  acquisition  for  p*  projection. 

n  Number  of  radians  in  a  semi-circle. 

(pctc  Reconstructed  phase  of  subject  motion. 

4>mt  Phase  of  scanner  rotation. 

^sub  Phase  of  subject  motion. 

I.  Introduction 

EVERAL  imaging  modalities  are  available  to  provide 
four-dimensional  descriptions  of  transient  tissue  defor¬ 
mations  (i.e.,  three  spatial  dimensions,  one  temporal  dimen¬ 
sion)  [1].  X-ray  computed  tomography  (CT)  is  often  used  in 
this  regard  due  to  its  excellent  spatial  resolution,  despite  the 
risks  associated  with  the  use  of  ionizing  radiation  [2].  However 
scanner  rotation  speed  limits  temporal  resolution  using  stan¬ 
dard  image  reconstruction  algorithms  [3].  This  may  limit  the 
quality  of  images  during  dynamic  subject  motion,  especially 
for  deformations  at  high  frequencies.  As  a  result,  temporally- 
resolved  imaging  during  rapid,  periodic  motion,  e.g.,  high- 
frequency  and  multi-frequency  oscillatory  ventilation  [4],  [5], 
is  difficult  to  acquire.  Especially  with  faster  rotation  times, 
larger  axial  coverage,  and  the  ability  to  significantly  reduce 
radiation  exposure  utilizing  newly  emerging  CT  technolo¬ 
gies  [6],  [7],  there  is  a  renewed  interest  in  assessing  dynamic 
physiologic  phenomena  via  CT. 

With  CT  imaging,  the  X-ray  source  and  detector  array  rotate 
about  the  subject,  while  X-rays  are  projected  at  varying  angles. 
The  resulting  projections  on  the  opposing  detector  array  are 
characterized  by  linear  attenuation  coefficients,  representing 
the  tendency  of  various  biological  tissues  to  absorb  or  scatter 
X-rays  [8].  The  raw  data  acquired  by  the  scanner  comprises 
projections  at  multiple  angles.  A  two-dimensional  image  of  the 
spatially-distributed  linear  attenuation  coefficients  can  then  be 
obtained  by  reconstruction  of  the  projection  data. 

Various  reconstruction  algorithms  are  available  for  clin¬ 
ical  purposes,  such  as  filtered  back-projection  or  iterative 
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reconstruction  [9],  with  each  achieving  a  particular  balance 
between  image  quality  and  computational  cost.  Standard 
reconstruction  algorithms  assume  that  the  subject  is  stationary 
during  scanner  rotation,  such  that  each  acquired  projection 
captures  the  same  spatial  orientation  through  a  different  angle. 
However  if  the  subject  is  moving  during  image  acquisition, 
the  projection  data  will  contain  multiple  spatial  orientations 
and  projection  angles  [10].  In  such  situations,  standard  recon¬ 
struction  algorithms  will  blur  the  resulting  image  due  to 
motion  artifact  [11].  This  becomes  especially  problematic  dur¬ 
ing  repetitive  periodic  motion  of  anatomic  structures  at  rates 
equal  to  or  faster  than  the  scanner  rotation  frequency  [4],  [12]. 
Motion  artifact  mitigation  and/or  temporally  resolved  image 
reconstruction  in  three  dimensions  is  possible  using  appro¬ 
priate  gating  or  triggering  techniques  [12].  However  such 
techniques  selectively  target  image  acquisition  at  the  same 
subject  state  repeatedly  over  the  course  of  multiple  motion 
cycles,  rather  than  acquiring  multiple  transient  states  within 
one  cycle. 

Image  acquisition  during  periodic  subject  motion  has  been 
well-described  in  the  fields  of  cardiac  and  thoracic  imaging, 
for  which  cardiogenic  motion  (in  the  frequency  range  of 
1  to  2  Hz)  and  spontaneous  breathing  (0.1  to  0.7  Hz)  present 
sources  of  motion  artifact.  Gating  techniques  for  such  appli¬ 
cations  may  utilize  external  signals  for  either  retrospective 
gating  or  prospective  triggering  of  the  image  acquisition 
cycle.  In  thoracic  CT,  a  variety  of  external  signals  including 
respiratory  flow,  pressure,  or  abdominal  motion  may  be  used 
to  target  image  timing  according  to  lung  volume  or  respira¬ 
tory  phase  (i.e.,  end-expiration  and  end-inspiration)  [13]-[15]. 
In  cardiac  CT,  an  electrocardiogram  (EKG)  may  be  used  to 
target  image  timing  according  to  cardiac  phase  and  heart 
rate,  identifying  a  period  for  which  the  mediastinal  contents 
remain  relatively  stationary  (i.e.,  diastole)  [12],  [16]-[18]. 
However  conditions  involving  highly-dynamic  subject  states, 
such  as  high-frequency  oscillatory  ventilation  (HFOV),  do  not 
provide  comparable  motionless  time  windows.  HFOV  is  an 
alternative  form  of  mechanical  ventilation  that  may  be  used 
during  severe  refractory  hypoxemia  in  the  Acute  Respira¬ 
tory  Distress  Syndrome  (ARDS).  In  contrast  to  conventional 
forms  of  mechanical  ventilation,  HFOV  delivers  tidal  volumes 
that  may  be  smaller  than  anatomic  dead  space,  at  rates 
of  2  to  20  Hz  [19].  Thus  imaging  during  HFOV  has  been 
limited  to  planar  views  [20]  or  phase-averaged  volumet¬ 
ric  images  [4],  [21]-[23],  resulting  in  poor  spatio-temporal 
resolution. 

To  achieve  high-resolution  CT  images  during  such  highly 
dynamic  states,  an  alternative  approach  must  be  imple¬ 
mented  that  does  not  rely  on  the  reconstruction  of 
sequentially  acquired  projections  triggered  by  physiological 
events  [24],  [25].  In  this  study,  we  propose  the  use  of 
frequency-selective  4D-CT  (FSCT)  for  capturing  periodic  sub¬ 
ject  motion  with  both  high  spatial  and  temporal  resolutions. 
Our  approach  relies  on  the  continuous  acquisition  of  many 
projections  during  multiple  scanner  rotations  over  a  specified 
length  of  axial  coverage,  followed  by  a  retrospective  analysis 
of  the  resulting  phase  distribution  of  subject  motion  obtained 
from  projection  data. 


II.  Methods 

A.  Theory 

Consider  an  X-ray  source  rotating  through  a  circular  tra¬ 
jectory  at  frequency  /mt,  with  an  opposing  detector  array 
acquiring  projections  at  varying  angles,  during  periodic  subject 
motion  at  frequency  /sub-  The  phase  of  subject  motion  /sub 
can  be  expressed  as  a  function  of  time: 

/sub  (0  =  (2;r/sub0  mod  In,  t  >0  (1) 

where  each  distinct  value  of  /sub  between  0  and  'In  corre¬ 
sponds  to  a  particular  spatial  orientation,  cyclically  repeated 
during  the  periodic  motion.  Similarly,  the  phase  of  the  scanner 
rotation  is: 


/rot  (0  =  (27r/rotO  mod  2n,  t  >0  (2) 


where  /mt  is  equivalent  to  the  projection  angle  relative  to  the 
initial  projection  angle  at  time  f  =  0.  To  reconstruct  a  “still” 
image  of  the  subject  at  one  particular  phase  /rec,  projections 
must  be  acquired  at  times  tk  such  that  /sub  (tk)  =  /rec^ 


tk 


/rec  T  27r  k 
2^7  /sub 


where  A:  =  0,  1,  2, . . . 


(3) 


The  relative  X-ray  projection  angle  at  each  tk  is  given  by: 


where 


/rot  (tk)  — 


(/rec  +  2nk) 


/:  =  0,1,2,... 


mod  2n, 


(4) 


Multiple  distinct  “still”  images  of  the  subject  throughout  the 
periodic  motion  may  be  reconstructed  by  first  partitioning 
subject  phase  into  Vbin  equally  spaced  bins  on  the  interval 
0  to  2;r: 

1 


^  bin  ■ 


/Vbi, 


2n 

2n 

2n 


0 

<  / 

^bin 

1 

<  / 

-^bin 

^bin 

- 1 

/^bin 

2 

/Vbin 


Nbh 


^  /sub  <  2;r 


Nbh 


(5) 


and  then  grouping  the  sequentially  acquired  projections 
according  to  the  corresponding  subject  phase  bin: 


mod  Vbin  , 


f  >  0  (6) 


where  [  J  denotes  rounding  the  enclosed  argument  to  the 
largest  integer  less  than  or  equal  to  itself.  Equation  (6)  can  be 
rewritten  as: 


^bin  (t)  —  [(/sub^/Vbin)  mod  VbinJ  ,  t  >  0  (7) 

Binning  the  acquired  projections  according  to  subject  phase 
yields  Abin  sinograms,  each  containing  a  sufficient  number  of 
angles  required  for  image  reconstruction.  The  phase-binned 
sinograms  may  then  be  reconstructed  independently  to  produce 
Abin  images  in  a  sequence  representing  dynamic  subject  state 
during  periodic  motion.  Each  image  represents  subject  state 
over  a  limited  range  of  motion  characterized  by  the  size  of 
each  phase  bin.  Increasing  Abin  reduces  the  size  of  each 
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Fig.  1.  Illustration  of  FSCT  concept  for  arbitrarily  selected  parameters,  (a)  Scanner  phase  ^rot  as  a  function  of  time  for  a  2.0  FIz  rotation  frequency,  with 
25  projections  acquired  per  rotation  as  indicated  by  numbered  asterisks.  Projections  are  numbered  according  to  chronological  order  of  acquisition. 
Subject  phase  as  a  function  of  time  for  a  13.13  FIz  periodic  motion  cycle.  Subject  phase  is  labeled  according  to  3  equally  sized  bins,  indicated 
by  color  (red,  orange,  yellow),  (b)  Standard  CT  reconstruction  is  performed  on  sequentially  acquired  projections,  averaging  together  variable  subject 
phase  acquired  throughout  one  scanner  rotation,  (c)  FSCT  reconstruction  is  performed  on  phase-binned  subsets  of  projections,  each  containing  a 
limited  range  of  subject  phase. 


phase  bin,  thereby  reducing  motion  artifact  in  the  reconstructed 
images. 

Fig.  la  shows  an  example  of  transient  variations  in  subject 
motion  phase  during  scanner  rotation,  according  to  (1)  and  (2). 
Standard  CT  reconstruction  is  performed  on  sequentially 
acquired  projections,  resulting  in  variable  subject  motion  phase 
obtained  during  one  scanner  rotation,  illustrated  in  Fig.  lb. 
FSCT  reconstruction  is  instead  performed  on  subsets  of  pro¬ 
jection  data  binned  according  to  subject  motion  phase.  Fig.  Ic 
illustrates  more  uniform  subject  motion  phase  obtained  in  each 
binned  subset  of  projection  data  using  the  FSCT  approach 
described  in  (5). 

B.  Feasible  Sampling  Criterion 
In  practical  applications,  the  number  of  unique  projection 
angles  is  given  by  the  number  of  projections  acquired  per 
scanner  rotation  A^pro-  Projections  are  acquired  at  a  fixed 
rate  and  with  equal  angular  spacing  throughout  each  scanner 
rotation.  Therefore  projection  angle  (prot  may  be  considered 
a  discrete  variable  rather  than  a  continuous  function  of  time, 
such  that  (prot  for  the  p*  sequentially  acquired  projection  is 
given  by: 

2np 

— mod  277,  p  =  0,1,2,...  (8) 

-^pro 

Each  unique  value  of  cpmx  ip)  may  be  indexed  by  such 
that: 

«pro  ip)  =  P  mod  Afpro,  p  =  0,  1, 2, . . .  (9) 


Acquisition  of  the  p*  projection  occurs  at  time  point  tp: 

P  =  0,l,2,...  (10) 

^pro/rot 

with  subject  motion  phase: 

^sub  ip)  =  P  mod  In,  p  =  0,1,2,...  (11) 

-^pro/rot 

and  subject  phase  bin  index  /ibin^ 


«bin  (P)  = 


A^bin  /sub  ,  , , 
p -  mod  A'bi: 


(12) 


Each  projection  is  therefore  assigned  a  pair  of  indexes 
(npro,  «bin)  characterizing  the  associated  projection  angle  and 
subject  motion  phase  bin.  All  unique  combinations  of  npro  and 
nbin  must  be  sampled  to  appropriately  reconstruct  the  complete 
image  sequence  with  minimal  spatial  aliasing,  although  it  is 
possible  to  poorly  select  /sub,  /rot,  and  Abin  such  that  the 
acquired  projections  repeat  non-unique  pairs  (npro(p),  nbin(p)) 
before  complete  sampling.  The  relative  periodicity  between 
scanner  rotation  and  subject  motion  is  determined  by  the 
number  of  projection  acquisitions  which  may  occur  before  a 
repetition  in  the  pair  (/rot(p),  /sub(p)),  hereafter  referred  to  as 
Rpxo-  Projection  angle  index  npro  is  periodic  in  p  with  period 
Apro,  therefore  Rpm  must  be  an  integer  multiple  of  Npm'- 


Rpm  —  ^rotA^pro 


(13) 


where  R^x  is  the  smallest  integer  number  of  rotations  which 
may  complete  before  a  repetition  in  {/rot(p),  /sub(p))  occurs. 
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Fig.  2.  The  space  of  feasible  sampling,  showing  the  combinations  of  A/t,jn  and  frequency  ratio  (fgub/ W)  which  the  projection  angles  and  subject 
phase  bins  are  either  completely  (white)  or  incompletely  (black)  sampled  before  repetitious  sampling.  Parameter  selection  along  the  black  lines  may 
not  provide  sufficient  sampling  for  image  reconstruction.  Typical  ranges  of  fgub  ^tid  ffot  corresponding  to  conventional  mechanical  ventilation  (CMV), 
spontaneous  breathing,  cardiogenic  motion,  and  high-frequency  oscillatory  ventilation  (HFOV)  are  shown. 


Substituting  this  expression  for  Rp^  for  p  in  (11)  yields: 

4>sub  {Rpro)  =  mod  271  (14) 

/rot 

Thus  the  pair  {/'rot(f’), /’sub(p))  repeats  non-unique  combi¬ 
nations  when  an  integer  number  of  scanner  rotations  have 
completed  such  that: 

^sub  =  (15) 

/rot 

where  7?sub  is  the  smallest  integer  number  of  motion  cycles  that 
may  complete  before  a  repetition  in  {cprotip),  />sub(p))  occurs. 
If  /sub  and  /rot  are  incommensurable,  then  no  finite  values  of 
Rmt  and  Rsub  exist.  Instead,  both  R^t  and  R^ub  will  be  infinite, 
since  an  exact  repetition  in  (/rot(p),  /subCp))  will  never  occur. 
However  if  /sub  and  /rot  are  commensurable,  then  R^t  may 
be  found  by  continued  fraction  expansion  [26]. 

Finally  the  number  of  unique  pairs  (/rotCp),  /subCp))  for  a 
given  /rot  and  /sub  is  limited  by  an  upper  bound,  determined 
as  the  number  of  projections  acquired  per  scanner  rotation 
Afpro  multiplied  by  the  number  of  scanner  rotations  which 
may  occur  before  repetition  /?rot  ■  To  ensure  that  undersampling 
does  not  occur,  the  product  Afpro  ■  Rmt  must  be  greater  than 
the  total  number  of  samples  to  be  acquired,  which  is  given  by 
the  product  Np^  ■  Nun- 

Afnro  R^ot  f^rot 

P  =  — IlL  >  1  (16) 

iVproAfbin  M)m 

Note  that  this  feasible  sampling  criterion  does  not  depend  on 
Npm,  but  rather  only  on  Nbin,  /rot,  and  /sub-  Fig-  2  illustrates 
the  space  of  feasible  sampling,  given  a  selection  of  Nbin  and 
frequency  ratio  (/sub//rot)-  This  fractal-like  criterion  is  charac¬ 
terized  by  the  Thomae  function,  a  modified  Dirichlet  function 
that  is  non-zero  only  at  rational  values  of  its  argument. 

C.  Optimal  Sampling  Condition 

The  minimum  scanning  duration  for  a  feasible  selection  of 
Nbin  and  frequency  ratio  (/sub //rot)  depends  on  the  length 
of  time  required  to  completely  sample  all  unique  pairs  of 
phase  bin  and  projection  angle.  Under  ideal  circumstances. 


each  newly  acquired  projection  corresponds  to  a  previously 
unsampled  pair  (nproCp),  nbinCp)),  i-C-,  with  minimal  redun¬ 
dancy.  In  this  case,  the  scanning  duration  measured  in  number 
of  scanner  rotations  (Nmt)  is  lower-bounded  by  the  number  of 
motion  phase  bins: 


Nmt 

Nbin 


>  1 


(17) 


However  (17)  assumes  that  the  acquired  projections  uniformly 
sample  both  subject  phase  bins  and  projection  angle  bins 
without  redundancy.  More  generally,  the  degree  of  redun¬ 
dant  sampling  depends  on  the  frequency  ratio  (/sub//rot), 
as  indicated  by  Fig.  2.  Due  to  the  discrete  nature  of  phase- 
binning,  redundant  sampling  of  indices  {ripm^p),  HbinCp))  may 
occur  before  an  exact  repetition  of  phases  (npro(p),  «bin(7’)) 
is  encountered.  For  any  choice  of  Aj,in,  the  ideal  sampling 
condition  given  by  (17)  is  ensured  with  equality  when: 

^  —  (18) 

/rot  Nbin 

where  c  is  any  positive  integer  such  that  c  and  A),in  are 
coprime.  If  (18)  is  satisfied,  then  Nbin  scanner  rotations  may  be 
used  to  reconstruct  Nbin  motion  phase  bins.  Some  examples  of 
optimal  frequency  ratios  for  desired  Aj,in  values  are  provided 
in  Table  I,  according  to  (18).  Deviating  from  the  frequency 
ratios  given  by  (18)  may  result  in  substantially  increased 
scanning  duration  required  to  properly  reconstruct  the  image 
sequence. 


D.  Simulated  Phantom 

An  oscillating  imaging  phantom  was  simulated  using  a 
Shepp-Logan  design  [27],  with  two  ellipses  modified  by  sinu¬ 
soidal  oscillations  in  size,  intensity,  and/or  centroid  position 
at  a  rate  of  9.924  Hz,  as  shown  in  Fig.  3.  Image  size  was 
256  by  256  pixels.  An  X-ray  CT  scanner  was  simulated  with 
1000  detectors,  1600  projections  per  rotation,  and  3.509  Hz 
rotation  frequency  (285  ms  period).  These  frequencies  were 
chosen  to  maximize  the  number  of  feasible  choices  for 
Abin.  Projection  acquisition  was  performed  using  the  Radon 
transform  in  MATLAB  (The  Mathworks  Inc.,  Natick,  MA), 
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TABLE  I 

Optimal  Sampling  Conditions  to  Minimize  Scanning  Duration 


^bin 

Optimal  /sub //rot  ratio 

1 

Any 

2 

1/2  ,  3/2  ,  5/2  ,  7/2 , 9/2  , 

3 

1/3  ,  2/3 , 4/3  ,  5/3  ,  7/3  ,  ... 

4 

1/4 , 3/4 , 5/4 , 7/4 , 9/4  ,  . . . 

5 

1/5 , 2/5 , 3/5 , 4/5 , 6/5  ,  ... 

^bin 

c/lVbin  ,  with  c  and  IVbin  coprime 

Under  such  conditions,  required  scanning  duration  iV^ot  =  iV|,in. 


assuming  parallel  beam  projection  and  a  flat  detector  array. 
A  total  of  224,000  projections  over  140  scanner  rotations  were 
acquired  during  continuous  phantom  oscillation. 

Image  reconstruction  was  performed  by  Altered  back- 
projection  using  the  inverse  Radon  transform,  with  linear  inter¬ 
polation  and  Ram-Lak  ramp  Altering.  Standard  reconstruction 
was  performed  using  either  1  complete  rotation  (1600  sequen¬ 
tially  acquired  projections),  or  140  rotations  (224,000 
sequentially  acquired  projections).  FSCT  reconstruction  was 
performed  by  phase-binning  the  projection  data  according  to 
the  9.924  Hz  phantom  cycle,  using  a  range  of  Nmt  between 
1  and  140  rotations,  and  a  range  of  Abin  between  1  and 
100  motion  phase  bins.'  FSCT  reconstruction  error  was 
computed  as  the  mean  and  standard  deviation  of  voxel-wise 
absolute  intensity  difference  between  the  reconstructed  image 
and  the  phantom  image  across  all  reconstructed  phase  bins. 
The  reference  phantom  image  for  each  phase  bin  was  selected 
at  the  motion  phase  correspondingto  the  midpoint  of  the 
phase  bin. 

E.  Thoracic  Imaging  in  vivo 

One  pig  weighing  45  kg  was  scanned  with  approval  from  the 
University  of  Iowa  Institutional  Animal  Care  and  Use  Com¬ 
mittee  (Protocol  number  5031314).  The  animal  was  sedated 
using  an  intravenous  infusion  of  propofol  at  9  mg  kg“* 
hr“*,  and  paralyzed  with  intermittent  doses  of  rocuronium  at 
1  mg  kg“*.  EKG  waveforms  were  obtained  using  a  Philips 
patient  monitor  equipped  with  the  M3001A  measurement 
module  (Philips  Healthcare,  Andover,  MA).  The  animal  was 
ventilated  by  a  FabianHFO  ventilator/oscillator  (Acutronic 
Medical  Systems  AG,  Hirzel,  Switzerland)  using  pressure- 
controlled  ventilation  at  a  rate  of  24  min“'  with  5  cmH20 
positive  end-expiratory  pressure,  20  cmH20  peak  inspiratory 
pressure,  and  I:E  ratio  1:2.  The  animal  was  positioned  in 
the  gantry  of  a  SOMATOM  Eorce  CT  scanner  (Siemens 
Healthineers,  Eorchheim,  Germany)  such  that  the  axial  field  of 
view  was  centered  approximately  3  cm  caudal  to  the  carina. 
Scanning  was  performed  during  1)  ventilation  with  pressure- 
controlled  conventional  mechanical  ventilation  (CMV);  2) 
high-frequency  oscillatory  ventilation  (HFOV)  set  at  5  Hz; 
3)  HFOV  set  at  20  Hz;  and  4)  multi-frequency  oscillatory 
ventilation  (MFOV)  using  a  broadband  excitation  waveform 
containing  energy  set  at  5  Hz,  10  Hz,  15  Hz,  and  20  Hz  [5]. 
Each  of  the  oscillatory  waveforms  was  delivered  at  10  cmH20 

*A  sample  MATLAB  script  for  this  computational  phantom  simulation  is 
available  from  the  authors  upon  request. 


mean  airway  pressure,  100  cmH20  peak-to-peak  pressure, 
and  10  L  min“'  bias  flow.  The  animal  was  continuously 
scanned  at  80  kVp  tube  voltage  and  150  mA  tube  current 
for  a  total  duration  of  approximately  30  seconds,  at  a  4.0  Hz 
scanner  rotation  frequency  (250  ms  period).  The  Force  scanner 
acquired  1120  projections  per  rotation  using  920  detectors 
per  row  and  96  detector  rows.  Transverse  beam  width  was 
0.48  mm  at  the  isocenter,  with  0.6  mm  spacing  between 
detector  rows.  As  part  of  another  study  under  this  protocol, 
a  total  of  177  mL  Isovue-370  contrast  agent  (Bracco  Diagnos¬ 
tic  Inc.,  Monroe  Township,  NJ)  was  given  over  approximately 
one  hour.  Scanning  for  this  study  was  performed  20  minutes 
following  the  final  54  mL  dose. 

The  ventilatory  and  cardiac  frequencies  were  estimated 
from  the  ventilator  settings  and  EKG  heart  rate  monitor, 
however  these  estimates  lacked  sufficient  precision  for  the 
FSCT  reconstruction  technique.  Since  the  in  vivo  projection 
data  contained  substantial  frequency  content  at  both  ventilatory 
and  cardiac  frequencies,  temporal  variations  in  each  set  of 
projection  data  were  examined  to  determine  frequency  content 
prior  to  selection  of  FSCT  reconstruction  frequencies.  The 
acquired  projections  were  re-sampled  at  intervals  of  In  radians 
in  projection  angle,  producing  subsets  of  the  projection  data 
containing  transient  variations  in  subject  position  sampled 
at  the  frequency  of  the  scanner  rotation  (/rot)-  In  the  two- 
dimensional  discrete  Fourier  transform  (2D-DFT)  of  this  sub¬ 
set  of  projection  data,  the  dimensions  of  the  transformed 
projection  data  correspond  to  detector  position  and  time. 
Integrating  the  magnitude  of  the  2D-DFT  along  the  dimension 
corresponding  to  detector  position  produced  a  spectrum  of 
oscillatory  energy  in  the  projection  data  as  a  function  of 
temporal  frequency,  relative  to  scanner  rotation  frequency 
(i.e.,  frequencies  greater  than  0.5/rot  will  be  aliased).  Spectral 
peaks  identified  in  the  oscillatory  energy  spectrum  correspond 
to  either  ventilatory  or  cardiogenic  sources  of  motion.  Thus 
ventilatory  and  cardiac  frequencies  were  determined  from 
the  projection  data,  using  a  priori  knowledge  obtained  from 
physiologic  monitors  (e.g.,  ventilator  settings,  EKG-derived 
heart  rate)  to  determine  de-aliased  frequency  content. 

FSCT  reconstruction  was  performed  by  first  phase-binning 
the  projection  data  according  to  either  ventilatory  or  car¬ 
diac  cycle  as  described  in  the  Theory  section,  and  sub¬ 
sequently  reconstructing  each  phase-binned  sinogram  using 
a  Feldkamp-type  cone-beam  reconstruction  algorithm  [28] 
(Michigan  Image  Reconstruction  Toolbox,  J.  A.  Fessler,  Uni¬ 
versity  of  Michigan,  Ann  Arbor,  MI)  and  MATLAB  (The 
Mathworks  Inc.,  Natick,  MA). 

III.  Results 

Fig.  4  shows  the  simulated  phantom  at  time  f  =  0  seconds, 
alongside  standard  CT  reconstructions  and  FSCT  reconstruc¬ 
tion  of  projection  data  acquired  over  140  scanner  rotations. 
The  FSCT  image  is  only  shown  at  the  first  of  10  motion  phase 
bins.  A  transient  profile  is  shown  through  a  single  row  from 
the  phantom  and  FSCT  images,  resolving  the  periodic  motion 
of  the  phantom  over  multiple  9.924  Hz  oscillation  cycles. 
A  static  profile  through  the  same  row  is  provided  for  the 
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Fig.  3.  Simulated  phantom  shown  at  1 0  discrete  phases  equally  spaced 
throughout  a  continuous  and  periodic  deformation  cycle.  Two  ellipses  of 
the  standard  Shepp-Logan  phantom  [27]  were  modified  with  sinusoidal 
oscillations  in  size,  intensity,  and/or  centroid  position.  Image  size  was 
256  X  256  pixels. 

standard  reconstructed  images.  The  standard  reconstructions 
are  affected  by  motion  artifacts  that  obscure  hner  details  in 
lung  structure,  whereas  the  FSCT  reconstruction  yields  crisp 
spatial  boundaries  within  each  bin  and  smooth  transient  motion 
between  bins. 

Fig.  5  shows  the  error  of  FSCT  reconstruction  performed 
using  varying  subsets  of  the  total  acquired  projection  data  and 
varying  phase  bin  sizes.  Error  was  assessed  by  the  mean  and 
standard  deviation  of  voxel-wise  intensity  difference  between 
each  reconstructed  image  compared  to  the  reference  phantom 
images.  The  lower  bound  for  adequate  scanning  duration 
given  by  (17)  is  indicated  by  the  dashed  line,  delineating 
two  domains  with  different  behavior  of  the  reconstruction 
error  with  respect  to  the  amount  of  sampling  and  number  of 
bins.  The  minimum  sampling  criterion  is  not  satisfied  at  short 
scanning  durations  and  large  number  of  bins,  thus  aliasing  and 
streak  artifact  caused  by  undersampling  produce  substantial 
reconstruction  errors.  The  minimum  sampling  criterion  is 
satished  at  long  scanning  durations  and  small  number  of  bins, 
with  mean  voxel-wise  intensity  difference  converging  with 
increased  sampling.  However  when  A/bin  becomes  too  small, 
the  reconstruction  error  increases  as  a  result  of  the  increased 
bin  size  and  motion  artifact.  For  reference,  the  standard 
reconstruction  using  1  rotation,  shown  in  Fig.  4e,  corresponds 
to  the  location  N^in  =  1  and  A/rot  =  1  in  the  bottom-left  corner 
of  Figs.  5a  and  5b.  Likewise,  the  standard  reconstruction 
using  140  rotations,  shown  in  Fig.  4g,  corresponds  to  the 
location  A^bin  =  1  and  Nmt  =  140  in  the  top-left  corner 
of  Figs.  5a  and  5b. 

A.  Thoracic  Imaging  in  vivo 

Fig.  6  shows  frequency  content  in  the  projection  data 
acquired  during  in  vivo  mechanical  ventilation.  The  spectral 
peaks  correspond  to  either  cardiogenic  or  ventilatory  oscil¬ 
lations,  potentially  aliased  relative  to  the  scanner  rotation 
frequency.  FSCT  reconstruction  frequencies  were  chosen  for 
each  set  of  projection  data  based  on  the  spectral  peaks,  using 
a  priori  knowledge  to  de-alias  higher  frequency  content.  Fig.  7 
shows  example  images  reconstructed  at  both  cardiac  and 
ventilatory  frequencies  using  the  FSCT  approach.  Transient 
profiles  from  each  image  are  presented  through  a  single 
slice  from  the  transverse  plane,  resolving  the  motions  of 
various  intrathoracic  organs  and  tissues  across  subject  phase 


bins,  equally  spaced  throughout  the  periodic  motion  cycle. ^ 
Motion  artifact  is  present  in  each  image  wherever  extraneous 
motion  occurs  at  non-harmonic  frequencies  relative  to  the 
reconstruction  frequency — for  example,  cardiogenic  motion 
produces  blurring  in  images  reconstructed  at  the  ventilatory 
frequency,  and  vice  versa. 

It  should  be  emphasized  that  each  pair  of  ventilatory  and 
cardiac  reconstructions  presented  in  Fig.  7  were  obtained 
from  the  same  projection  data.  Furthermore  the  transient 
profiles  are  shown  through  the  same  spatially-oriented  slice. 
Therefore,  any  differences  in  image  quality  or  structural  clarity 
between  these  paired  reconstructions  are  a  direct  consequence 
of  modulating  the  reconstruction  frequency  and/or  number  of 
phase  bins. 

IV.  Discussion 

In  this  study,  we  have  presented  an  algorithm  for  recon¬ 
structing  temporally-resolved  CT  images  of  periodic  sub¬ 
ject  motion.  This  frequency-selective  4D-CT  (FSCT)  imaging 
technique  involves  continuously  acquiring  X-ray  projections 
during  multiple  rotations  of  the  source  around  the  oscillating 
subject,  then  sorting  the  acquired  projections  into  multiple 
subject  phase  bins,  each  yielding  a  complete  sinogram  for 
reconstruction  of  an  image  sequence. 

Respiratory-gated  and  non-gated  4D-CT  are  suitable 
for  use  during  irregular,  spontaneous  breathing  with  low 
subject  motion  frequency  compared  to  scanner  rotation 
frequency  [10],  [29].  FSCT  extends  4D-CT  imaging  across  a 
much  wider  range  of  application  in  terms  of  subject  motion 
frequencies,  yet  imposes  a  restricting  assumption  of  regular, 
periodic  motion.  Nevertheless,  the  limitations  imposed  by 
the  feasible  sampling  criterion  (16)  do  not  restrict  practical 
4D-CT  reconstruction  at  low  frequencies  such  as  spontaneous 
breathing.  This  is  indicated  in  Fig.  2  by  the  available  selection 
of  Nbin  up  to  20  for  low  /sub  relative  to  /rot-  However  at 
high  frequencies  of  subject  motion,  the  choice  of  Nbin  may 
be  infeasible  depending  on  (16). 

Traditional  retrospective  gating  for  cardiac  CT  can  be 
performed  by  continuously  acquiring  both  EKG  and  X-ray 
projections  as  the  subject  slowly  translates  through  the  scan¬ 
ner  [30].  The  projection  data  is  then  selectively  filtered  prior 
to  reconstruction,  retaining  projections  corresponding  to  the 
diastolic  phase  and  discarding  the  rest.  Despite  the  potential 
for  high  radiation  exposure  during  the  scan,  a  large  portion 
of  the  projection  data  may  actually  not  be  used  for  the  image 
reconstruction.  Alternatively  with  prospective  triggering,  EKG 
detection  of  diastole  initiates  the  scanner  to  emit  and  detect 
X-ray  beams  through  a  subset  of  projection  angles,  typically 
covering  180°  plus  half  the  angle  of  the  detector  array  [31]. 
The  subject  is  then  translated  to  scan  a  different  transverse 
plane  upon  the  next  diastolic  triggering  event.  This  type  of 
acquisition  is  informally  known  as  “step-and-shoot”  in  regard 
to  the  repeated  process  of  “stepping”  to  the  next  transverse 
plane  followed  by  “shooting”  the  X-ray  beam  for  acquisi¬ 
tion.  Compared  to  retrospective  gating,  prospective  triggering 

^Animations  of  subject  motion  are  provided  as  supplementary  materials, 
available  online  in  the  supporting  documents  /  multimedia  tab. 
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Fig.  4.  Comparison  between  standard  CT  reconstruction  (e,g)  and 
FSCT  reconstruction  (c)  of  the  simulated  phantom  oscillating  at  9.924  Hz 
(a).  A  transient  profile  shown  through  a  single  row  (white  dashed  line) 
is  provided  for  the  simulated  phantom  (b)  and  FSCT  reconstructed 
image  (d),  resolving  periodic  motion  over  multiple  9.924  Hz  oscillation 
cycles.  A  static  profile  through  the  same  row  is  provided  for  the  standard 
reconstructed  images  (f,h). 


reduces  the  amount  of  ionizing  radiation  required  for  cardiac 
imaging.  For  either  of  these  EKG-based  methods,  the  temporal 
resolution  of  the  image  is  determined  by  the  duration  of 
acquisition  within  each  cardiac  cycle  [32].  In  other  words, 
temporal  resolution  depends  on  how  much  cardiac  motion  is 
averaged  in  the  reconstruction:  faster  scanner  rotation  reduces 
the  duration  of  acquisition  within  each  cardiac  cycle,  and  thus 
results  in  less  motion  artifact. 

Multi-segment  reconstruction  is  another  effective  strategy 
for  improving  the  temporal  resolution  of  a  thoracic  image.  This 
approach  involves  reducing  the  duration  of  image  acquisition 
each  cardiac  cycle  by  acquiring  a  subset  of  projection  data  (i.e., 
180°  plus  half  of  the  detector  fan  angle)  over  multiple  cardiac 
cycles  [33].  Multi-segment  reconstruction  results  in  a  more 
consistent  state  of  the  mediastinum  during  each  acquisition 
compared  to  standard  EKG-triggering  approaches.  However 
the  extended  scanning  protocol  increases  susceptibility  to 
image  corruption  from  other  extraneous  motions  such  as  respi¬ 
ration  [34]-[36].  Several  other  techniques  have  been  developed 


Fig.  5.  Contours  of  error  in  FSCT  reconstructed  image  of  simulated 
phantom  for  varying  selection  of  scanning  duration  and  number  of 
bins  A/tjjn.  Error  was  assessed  by  (a)  mean  and  (b)  standard  deviation 
of  voxel-wise  absolute  intensity  difference  between  the  reconstructed 
image  and  the  phantom  image  across  all  reconstructed  phase  bins.  The 
lower  bound  for  adequate  sampling  given  by  (17),  is  indicated  by  the 
dashed  line,  using  A/(jet  =  1000,  A/pro  =  1600,  ffot  =  3.509  Hz,  fgup  = 
9.924  Hz. 


for  mitigating  motion  artifact  and  radiation  exposure  in  cardiac 
imaging,  such  as  compensation  algorithms  that  modify  the 
projection  data  to  approximate  a  motionless  image  of  the 
thorax  based  on  estimated  motion  models  [37]-[41]. 

Alternative  strategies  for  volumetric  imaging  during 
dynamic  motion  using  X-ray  CT,  such  as  quasi¬ 
static  [42]  or  temporally-averaged  [4],  [43]  approaches, 
tolerate  some  amount  of  either  temporal  or  spatial  blurring 
when  using  sequential  projection  acquisition.  ESCT  fully 
resolves  the  spatio-temporal  dynamics  of  subject  motion,  yet 
only  throughout  the  cycle  of  one  specified  reconstruction 
frequency.  Thus,  ESCT  images  represent  a  periodic, 
dynamic  steady  state  of  the  subject.  Imaging  high-frequency 
motions  (i.e.,  greater  than  the  scanner  rotation  frequency)  is 
an  important  potential  application  of  the  ESCT  technique, 
opening  a  relatively  unexplored  area  of  research  in  CT 
imaging.  Reconstruction  algorithms  incorporating  direct 
sinogram  alterations  introduce  some  uncertainty  in  the 
resulting  image  due  to  the  use  of  numerically  estimated 
motion  models  used  to  approximated  the  actual  subject 
motion  [37]-[41].  Motion  estimation  approaches  maintain 
reasonably  low  exposure  to  ionizing  radiation.  Our  ESCT 
technique  is  intended  to  provide  high-fidelity  images  of  the 
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Fig.  6.  Frequency  content  in  the  projection  data  acquired  during  in  vivo 
mechanical  ventilation  of  a  pig,  using  either  conventional  mechanical  ven¬ 
tilation  (a),  high-frequency  oscillatory  ventilation  (b,c),  or  multi-frequency 
oscillatory  ventilation  (d).  Frequency  on  the  horizontal  axis  obtained 
from  discrete  Fourier  transform  using  the  scanner  rotation  frequency 
f|,Q(  =  4  FIz  as  the  sampling  frequency.  Thus  subject  motion  frequencies 
greater  than  O.Sfpot  are  aliased.  Spectral  peaks  correspond  to  either 
ventilatory  motion  or  cardiogenic  motion.  Spectral  peaks  corresponding 
to  heart  rate  are  indicated  by  asterisks. 


subject  during  actual  periodic  motion,  albeit  at  the  expense 
of  increased  radiation  exposure.  Thus  at  present  FSCT  may 
be  more  suitable  for  research  purposes,  rather  than  clinical 
diagnostic  use. 

The  selection  requirements  for  FSCT  imaging  parameters 
were  characterized  for  a  given  subject  motion  frequency 
/sub  and  desired  number  of  images  A^bin,  yielding  a  feasible 
sampling  condition  (16)  as  well  as  an  optimal  sampling 
condition  (18).  The  use  of  a  circular  X-ray  source  trajectory 
allows  improved  temporal  resolution  compared  to  similar 
phase-correlated  reconstruction  of  projections  acquired  during 
a  spiral  or  helical  trajectory  [24],  [25].  FSCT  imaging  is 
feasible  for  periodic  subject  motions  at  nearly  any  frequency 
provided  an  appropriate  frequency  ratio  (/sub//rot)  is  used  to 
ensure  full  sampling  without  repetition  (Fig.  2).  Thus  FSCT 
can  be  used  to  obtain  temporally-resolved  image  sequences 
of  both  high-  and  low-frequency  anatomic  motion.  Fig.  7 
demonstrates  several  successful  reconstructions  using  FSCT 
across  a  wide  range  of  frequencies  (0.4  FIz  to  20  Hz)  relative 
to  the  scanner  rotation  frequency  (4  Hz).  Ideally,  scanner 
parameters  such  as  fmt  could  be  tuned  to  minimize  feasible 
scanning  duration  given  a  specific  subject  motion  frequency 
/sub-  In  practical  applications,  most  commercially  available 
CT  scanners  offer  a  small  number  of  choices  for  fmt  and 
A^pro  preset  by  the  manufacturer.  However  unlike  conventional 
CT  imaging,  the  temporal  resolution  for  FSCT  is  determined 
not  only  by  the  rotation  speed  of  the  scanner,  but  also  by  a 
combination  of  imaging  and  reconstruction  parameters. 


Other  potential  uses  of  FSCT  span  a  variety  of  physiologic 
applications  involving  periodic  motion,  including  analysis 
of  airway,  lung  tissue,  or  blood  vessel  deformation,  tumor 
localization  and  tracking,  cardiac  wall  motion  abnormalities, 
as  well  as  joint  motion.  Future  work  may  also  involve  four¬ 
dimensional  image  registration  to  quantify  regional  strains  of 
thoracic  structure  during  periodic  motion  [42],  [44].  Four¬ 
dimensional  registration  approaches  identify  smooth  trans¬ 
formations  between  volumetric  images  through  a  tempo¬ 
ral  dimension  (or  in  this  case  specifically,  subject  motion 
phase)  [45],  [46].  A  periodic  constraint  may  be  applied  to 
FSCT-registration,  ensuring  continuity  and  smoothness  of  the 
computed  transformation  throughout  the  entire  reconstructed 
sequence,  including  the  wrap-around  from  the  final  image  to 
the  initial  image. 

A.  Limitations 

Some  limitations  of  our  technique  should  be  noted.  First,  the 
extensive  scanning  duration  required  for  FSCT  reconstruction 
exposes  the  subject  to  prohibitive  radiation  dosage,  i.e.,  up  to 
two  orders-of-magnitude  greater  than  standard  X-ray  CT  imag¬ 
ing.  Accordingly,  clinical  implementation  of  FSCT  imaging  in 
patients  may  not  be  warranted  without  substantial  reductions  in 
radiation  exposure.  Exact  control  of  scanner  rotation  frequency 
/rot  according  to  (18)  may  allow  minimization  of  redundant 
projection  acquisition  and  thus  minimization  of  scanning  dura¬ 
tion.  Iterative  reconstruction  of  each  phase-binned  sinogram 
may  allow  dose  reductions  at  the  expense  of  increased  com¬ 
putational  cost  [47].  It  may  be  possible  to  reduce  radiation 
dose  further  by  minimizing  the  number  of  acquired  projection 
angles  per  subject  phase  bin  using  half-scan  approaches  for 
cone-beam  reconstruction  [48],  [49].  Reductions  in  scanning 
duration  may  be  achievable  using  specialized  dual-source  CT 
scanners  equipped  with  two  independent  X-ray  sources  and 
corresponding  detector  arrays  positioned  at  different  projection 
angles  through  the  subject  [50],  [51].  Finally,  prospective  trig¬ 
gering  approaches  generally  reduce  radiation  exposure  com¬ 
pared  to  retrospective  gating  approaches  [52],  [53].  Prospective 
FSCT  imaging  may  also  be  possible  by  synchronizing  X-ray 
beam  actuation  to  specified  phases  of  a  particular  periodic 
subject  motion  cycle,  with  frequency  known  a  priori.  In  this 
way,  prospective  triggering  could  reduce  the  total  radiation 
dose  by  selectively  avoiding  radiation  exposure  during  non- 
reconstructed  phases. 

Although  the  FSCT  technique  was  evaluated  quantitatively 
in  computational  simulations,  evaluation  of  the  experimental! 
obtained  images  presented  in  Fig.  7  was  qualitative.  Quantita¬ 
tive  assessment  of  anatomic  structures  in  FSCT  images,  such 
as  segmented  volumes  of  intrathoracic  structures,  may  require 
“step-and-shoot”  implementation. 

Finally,  the  FSCT  technique  proposed  herein  uses  only  one 
value  of  /sub  per  image  sequence.  It  is  possible  to  reconstruct 
multiple  image  sequences  from  the  same  raw  projection  data, 
simply  by  using  a  different  value  of  /sub  for  each  sequence. 
Fig.  7  shows  several  cases  in  which  two  image  sequences  were 
reconstructed:  one  at  the  ventilatory  frequency  and  another 
at  the  cardiac  frequency.  However  the  resulting  images  will 
contain  resolved  motion  at  only  the  one  specific  reconstruction 
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Fig.  7.  Example  images  reconstructed  from  projection  data  acquired  over  30  seconds  of  scanning  during  in  vivo  mechanical  ventilation  of  a  pig, 
using  (a)  conventional  mechanical  ventilation  at  0.4  Hz  (24  min’^),  (b)  high-frequency  oscillatory  ventilation  at  5  Hz  (300  min’^),  (c)  high-frequency 
oscillatory  ventilation  at  20  Hz  (1 200  min’^ ),  and  (d)  multi-frequency  oscillatory  ventilation  using  a  broadband  waveform  containing  energy  at  5,  1 0, 
15,  and  20  Hz.  Both  ventilatory  and  cardiac  frequencies  were  reconstructed  from  the  same  projection  data.  Transient  profiles  from  each  image  are 
presented  through  a  single  slice  from  the  transverse  plane  over  1  to  6  seconds,  showing  resolved  motion  across  subject  phase  bins  equally  spaced 
throughout  the  periodic  motion  cycles. 


frequency,  as  well  as  any  harmonics  or  integer  multiples 
of  the  reconstruction  frequency.  All  other  extraneous  motion 
in  the  subject  will  cause  motion  artifact  and  blurring.  For 
this  reason,  blurring  around  the  heart  is  apparent  in  images 


reconstructed  at  the  ventilatory  frequency,  whereas  blurring  in 
the  lungs  is  apparent  in  images  reconstructed  at  the  cardiac  fre¬ 
quency.  It  may  be  possible  to  combine  information  from  both 
reconstructed  images  to  reduce  apparent  motion  artifact  [54]. 
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Additionally,  FSCT  may  not  be  appropriate  for  temporally 
resolving  physiologic  processes  that  are  either  irregular,  aperi¬ 
odic,  or  otherwise  occur  on  multiple  time  scales,  such  as  dif¬ 
fusive  equilibration  during  angiography  or  tracer  gas  washout. 

V.  Conclusion 

Frequency-selective  CT  imaging  produces  a  temporally- 
resolved  volumetric  image  sequence  of  periodic  subject  motion 
at  a  specified  reconstruction  frequency.  This  novel  technique 
enables  imaging  of  both  low-  and  high-frequency  dynamic 
periodic  motion  at  a  specified  frequency,  with  minimal  extra¬ 
neous  motion  artifact.  The  FSCT  approach  may  offer  improved 
imaging  fidelity  of  dynamic  subject  motion  compared  to  more 
conventional  ‘quasi-static’  approaches.  When  combined  with 
registration  techniques  [45],  FSCT  may  provide  detailed  four¬ 
dimensional  information  of  distributed  tissue  deformation  pre¬ 
viously  inaccessible  using  standard  reconstruction  approaches. 
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Introduction: 

Gas  flow  in  the  lung  during  oscillatory  ventilation  is  distributed  in  a  heterogeneous  and  frequency-dependent 
manner  [1],  However  severe  ventilation  heterogeneity  may  contribute  to  impairments  in  gas  exchange,  ventilation- 
to-perfusion  mismatch,  and  ventilator-induced  lung  injury  (VILI).  Minimizing  acinar  strain  heterogeneity  may 
reduce  the  risk  of  VILI  while  producing  efficient  carbon  dioxide  elimination.  In  this  study  we  used  a  gas  transport 
model  to  optimize  oscillatory  ventilator  waveforms,  in  which  the  spectral  content  of  each  flow  waveform  was 
adjusted  to  minimize  parenchymal  strain  heterogeneity  during  eucapnic  ventilation  [2,3]. 

Methods: 

A  heterogeneous  canine  lung  model  consisting  of  N  terminal  viscoelastic  acini  was  ventilated  with  a  simulated 
oscillatory  flow  waveform  composed  of  M  simultaneous  frequencies.  The  relative  magnitudes  of  flow  at  each 
frequency  was  numerically  optimized  according  to  a  cost  function  O,  defined  as  the  coefficient  of  variation  of 
acinar  volumetric  distension.  Using  a  Monte  Carlo  technique,  we  estimated  the  optimal  oscillatory  waveform  that 
produced  the  lowest  O,  and  thus  the  lowest  degree  of  parenchymal  strain  heterogeneity.  The  relative  magnitudes 
of  each  oscillatory  component  was  adjusted  to  achieve  eucapnia. 

Results: 

Optimal  oscillatory  waveforms  were  characterized  by  flows  with  large  amplitudes  at  low  frequencies,  combined 
with  small  amplitude  flows  at  high  frequencies.  Average  acinar  strain  during  eucapnic  ventilation  was  reduced 
when  additional  higher  frequency  components  were  included  in  the  waveform. 

Conclusion: 

Superposition  of  multiple  simultaneous  oscillatory  frequencies  provides  more  uniform  ventilation  distribution 
compared  to  single  frequency  oscillatory  ventilation,  as  well  as  more  mechanically  efficient  gas  exchange  to 
achieve  eucapnia.  An  optimal  combination  of  frequencies,  amplitudes,  and  phases  in  an  oscillatory  ventilator 
waveform  may  be  determined  according  to  the  distribution  of  flows  throughout  the  heterogeneous  lung  periphery. 
Work  supported  by  DOD  grant  PR151761 . 
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Introduction:  Oscillatory  flows,  when  used  in  a  lung  protective  ventilation  strategy,  are  distributed  throughout  the 
respiratory  system  in  a  heterogeneous  and  frequeney-dependent  manner  [1].  Ventilation  heterogeneity  contributes 
to  impairments  in  gas  exchange  and  ventilator-indueed  lung  injury  (VILI),  due  to  regional  overdistention  and 
atelectrauma.  Thus  minimizing  acinar  strain  heterogeneity  may  also  minimize  the  risk  of  VILI.  In  this  study,  we 
developed  an  optimization  algorithm  for  tuning  the  speetral  content  of  oseillatory  ventilator  waveforms,  in  order  to 
minimize  aeinar  strain  heterogeneity  in  eomputational  model.  Emphasis  was  placed  on  the  development  of 
oscillatory  waveforms  comprised  of  multiple  simultaneous  frequencies,  which  may  be  ‘tuned’  to  local  mechanical 
properties  [2]. 

Materials  &  Methods: 

For  optimization  of  M  simultaneous  frequencies  delivered  to  a  heterogeneous  lung  model  with  N  terminal 
viscoelastic  acini  [1],  we  defined  an  input  oscillatory  waveform  as  a  sum  of  sinusoids  with  M  discrete  frequencies 
and  amplitudes.  We  further  defined  an  objeetive  function  <1>,  as  the  eoeffieient  of  variation  of  acinar  root-mean- 
square  volumes: 


where  Vj-ms.n  is  the  root-mean-square  volume  over  time  in  the  n*  acinus,  and  is  the  average  across  all  N  acini. 
Minimization  of  O  therefore  minimizes  aeinar  strain  heterogeneity.  For  each  value  of  M  (i.e.  1,  2,  3,  ...),  100 
randomized  initial  guesses  based  on  a  Monte  Carlo  technique  were  used  for  a  numerieal  eonstrained  nonlinear 
minimization  routine  (MATLAB  v8.4.  The  Mathworks  Inc.,  Natick  MA). 

Results  &  Discussion:  The  optimal  solution  for  small  numbers  of  frequencies  (i.e.  M  <  4)  was  easily  determined 
using  a  random  search  over  the  design  variable  spaee.  Figure  1  shows  the  optimal  solution  found  using  random 
seareh  over  3  frequeneies.  The  optimal  solution  for  M  >  4  frequeneies  became  a  poorly  defined  problem  due  to 
the  presenee  of  many  loeal  minima  in  <1>.  Nonetheless,  the  optimal  solution  for  M  -I-  1  frequeneies  generally 
ineluded  the  same  frequeneies  found  in 
the  unique  solution  for  M  frequeneies, 
with  a  similar  amplitude  distribution. 

Conclusions:  Superposition  of  multiple 
simultaneous  frequencies  provides  more 
uniform  ventilation  distribution  compared 
to  single  frequeney  oseillatory  ventilation, 
with  less  potential  for  VILI.  An  optimal 
eombination  of  frequeneies,  amplitudes, 
and  phases  in  an  oscillatory  ventilator 
waveform  may  be  determined  according 
to  the  frequeney-dependent  distribution  of 
flows  throughout  a  heterogeneous  lung 
periphery. 

References: 
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Figure  1.  Convergence  of  the  minimization  algorithm  from  randomized 
initial  guesses  (indicated  by  grey  lines)  to  local  minima  (colored  spheres).  In 
this  case,  the  optimal  solution  (red  sphere)  forM  =  3  frequencies  (/i,  /2,  /s 
shown  in  Hz)  yielded  27%  reduction  in  O  compared  to  the  optimal  solution 
for  M  —  1  frequency.  The  resulting  optimized  flow  waveform  over  one 
second  is  shown. 
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Introduction:  Mechanical  ventilation  is  a  life-sustaining  intervention  delivering  flow  and  pressure  to  maintain  gas  exchange  in 
patients  with  acute  respiratory  failure.  Regional  maldistribution  of  ventilation  may  unintentionally  result  In  ventilatlon-to- 
perfusion  mismatching  or  ventilator-induced  lung  injury  (VILI),  especially  in  mechanically-heterogeneous  injured  lungs.  High- 
frequency  oscillatory  ventilation  (HFOV)  is  an  alternative  ventilation  modality  designed  to  minimize  VILI  through  the  delivery  of 
small  tidal  volumes  at  high  instantaneous  flow  rates.  However,  the  distribution  of  flow  throughout  the  lung  during  HFOV  may  be 
heterogeneous  and  frequency-dependent.  Multi-frequency  oscillatory  ventilation  (MFOV)  has  been  recently  proposed  to  enhance 
lung-protective  ventilation  by  delivering  flow  at  multiple  frequencies  simultaneously  [1].  In  this  study,  we  used  Frequency- 
Selective  CT  (FSCT)  imaging  [2]  combined  with  four-dimensional  (4D)  Image  registration  [3]  to  compare  the  regional 
distributions  of  periodic  structural  deformation  during  HFOV  and  MFOV. 

Materials  &  Methods:  Experiments  were  approved  by  the  University  of  Iowa  Institutional  Animal  Care  and  Use  Committee 
(Protocol  number  5031314).  Sedated  and  paralyzed  pigs  were  ventilated  using  conventional  mechanical  ventilation  (CMV)  at  a 
rate  of  20  min’*  and  tidal  volume  between  10-12  mL  kg’*,  HFOV  delivered  at  5  Hz,  and  MFOV  delivered  using  uniform  flow 
amplitudes  at  5,  10,  15,  and  20  Hz,  and  scanned  using  a  Siemens  Somatom  Force.  FSCT  was  used  to  reconstruct  temporally- 
resolved  Image  sequences  of  structural  deformations  during  periodic  ventilation  [2].  A  4D  tissue-volume  preserving  image 
registration  technique  was  then  used  to  obtain  regional  maps  of  expansion  and  contraction  relative  to  a  single  arbitrarily-selected 
reference  Image  [3].  Within  each  region-of-interest  in  the  reference  image,  periodic  expansion  and  contraction  were  assessed  by 
the  discrete  Fourier  transform  of  the  normalized  time-varying  Jacobian  determinants. 

Results  &  Discussion:  Total  volumetric  strain  (peak-to-peak  mean  +  standard  deviation)  during  HFOV  (7.9%  ±  3.1%)  and 
MFOV  (6.9%  +  2.5%)  were  substantially  lower  than  during  CMV  (31.3%  +/-  11.6%).  The  distribution  of  strain  amplitudes 
during  HFOV  delivered  at  5  Hz  tended  to  vary  primarily  In  the  ventral-dorsal  direction,  in  accordance  with  the  gravitational  field 
in  supine  position  (Figure  1).  Regional  strain  amplitudes  during  MFOV  varied  with  frequency,  with  different  lung  regions 
selectively  filtering  the  harmonic  frequency  content  of  the  broadband  oscillatory  flow. 
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Figure  1:  Regional  distributions  of  volumetric  strain  amplitudes  during  HFOV  and  MFOV  shown  in  transverse  cross- 
section,  using  discrete  Fourier  transform  of  normalized  time-varying  Jacobian  determinants. 

Conclusion:  Parenchymal  stretch  during  oscillatory  ventilation  is  regionally  heterogeneous  and  frequency-dependent.  The 
broadband  spectral  content  of  multi-frequency  oscillatory  ventilation  may  enhance  gas  transport  in  the  presence  of  periodic 
parenchymal  deformation  with  higher  harmonics,  which  may  be  further  adjusted  to  compensate  for  patient-specific  regional 
heterogeneity. 
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Background 

Primary  blast  lung  injury  (PBLI)  is  a  potentially  fatal  condition  affecting  victims  in  proximity  to 
explosive  detonations  in  combat  and  terrorist  attack  settings.  Sudden  stresses  and  strains 
throughout  the  lung  are  induced  by  pressure  waves,  resulting  in  parenchymal  tissue  damage  and 
derecruitment,  alveolar  hemorrhage,  air  emboli,  and/or  pneumothoraces.  The  resulting 
structural  injuries  may  lead  to  severe  ventilation-perfusion  mismatch,  hypoxia,  and  respiratory 
failure.  Although  mechanical  ventilation  is  often  required  in  blast  victims,  maldistribution  of  gas 
flow  in  a  mechanically  heterogeneous  lung  may  lead  to  further  injury.  In  this  study,  we  developed 
a  computational  model  of  a  human  lung  to  simulate  the  effects  of  PBLI  on  the  distribution  of  gas 
flow  during  mechanical  ventilation. 

Methods 

A  computational  model  of  a  human  airway  network  was  generated  using  central  airways 
obtained  from  segmented  X-ray  computed  tomographic  scans,  and  algorithmically  generated 
peripheral  airways  (M.H.  Tawhai  et  al.,  J  AppI  Physiol  97(6):2310-2321,  2004).  The  model 
consisted  of  60,494  cylindrical  airway  segments,  with  30,243  terminal  bronchi  subtended  by 
viscoelastic  acini.  Distributed  mechanical  properties  of  tissues  and  airways  were  simulated  to 
represent  typical  healthy  and  blast-injured  lungs,  the  latter  characterized  by  bilateral 
derecruitment  and  increased  tissue  stiffness,  focusing  on  the  perihilar  regions  of  the  model.  Gas 
flow  delivered  at  the  airway  opening  was  distributed  throughout  the  lung  according  to  regional 
mechanical  impedance,  at  frequencies  ranging  from  conventional  mechanical  ventilation  (0.2  Hz) 
to  high-frequency  oscillatory  ventilation  (2  to  20  Hz). 

Results 

Delivered  flow  was  heterogeneously  distributed  in  the  blast-injured  lung  during  both 
conventional  mechanical  ventilation  and  high-frequency  oscillatory  ventilation.  During 
conventional  mechanical  ventilation,  flow  was  distributed  primarily  according  to  local  tissue 
stiffness.  At  higher  frequencies,  the  distribution  of  flow  became  increasingly  heterogeneous  and 
frequency-dependent,  with  some  regions  being  underventilated  while  other  regions  experienced 
substantially  greater  distension. 

Conclusion 

The  distribution  of  parenchymal  distension  during  mechanical  ventilation  depends  on  injury 
severity  and  ventilatory  modality.  Computational  modeling  may  aid  in  predicting  the  progression 
of  ventilator-induced  lung  injury,  and  may  also  provide  a  means  for  investigating  alternative 
ventilation  modalities  employing  lung-protective  strategies. 
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Introduction:  Ventilator-induced  lung  injury  (VILI)  in  mechanically  ventilated  patients  may  result  from  extreme  stresses,  strains, 
and  strain  rates  associated  with  parenchymal  overdistension,  as  well  as  repetitive  alveolar  collapse  and  reopening.  VILI  is 
especially  hazardous  within  injured  lungs,  in  which  mechanical  heterogeneity  leads  to  heterogeneous  and  frequency-dependent 
distribution  of  peripheral  flows  and  pressures.  Multi-frequency  oscillatory  ventilation  (MFOV)  is  a  novel  lung-protective 
ventilatory  modality,  relying  on  the  delivery  of  oscillatory  flows  at  multiple  simultaneous  frequencies.  With  appropriate  spectral 
energy  content,  MFOV  waveforms  can  provide  more  even  ventilation  distribution  compared  to  traditional  single-frequency 
oscillatory  ventilation  (SFOV)  [1].  Since  the  goal  of  MFOV  is  to  improve  morbidity  and  mortality  in  ventilated  patients,  we 
hypothesized  that  the  spectral  content  of  MFOV  waveforms  could  be  tuned  to  minimize  risk  for  VILI.  To  test  this  hypothesis,  we 
developed  a  computational  lung  model  to  optimize  the  spectral  content  of  simulated  MFOV  waveforms.  Various  objective 
functions  were  used  to  modulate  the  relative  contributions  of  strain,  strain  rate,  and  mechanical  power  imparted  to  the 
parencyhmal  tissues  [2] . 

Methods:  The  computational  model  consisted  of  a  central  airway  tree  segmented  from  a  computed  tomographic  image  of  a 
supine  pig,  with  algorithmically  generated  peripheral  airways.  The  model  consisted  of  60,494  airway  segments  in  total,  with 
30,243  viscoelastic  acini.  Eucapnic  ventilation  was  simulated  using  a  Monte  Carlo  technique,  to  choose  from  MFOV  waveforms 
consisting  of  four  frequencies  (5,  10,  15,  and  20  Hz)  with  randomized  volume  amplitudes  and  phases.  Corresponding 
distributions  of  acinar  strain,  strain  rate,  and  mechanical  power  throughout  the  model  were  determined  for  each  MFOV 
waveform,  and  used  to  compute  objective  functions  with  variable  emphasis  on  the  relative  percent  contribution  to  VILI. 

Results:  Our  simulations  demonstrated  that  MFOV  waveforms  were  superior  to  SFOV  for  both  healthy  and  injured  lung  models. 
The  optimized  volume  amplitudes  at  each  frequency  in  the  MFOV  waveform  were  modulated  according  to  the  relative  objective 
function  weighting  of  strain  vs.  strain  rate,  as  well  as  the  degree  of  lung  mechanical  heterogeneity.  Increasing  contribution  of 
strain  rate  to  VILI,  as  opposed  to  strain,  resulted  in  decreasing  volume  amplitudes  at  higher  frequencies  in  the  optimized  MFOV 
waveforms.  Alternatively,  optimized  MFOV  waveforms  minimizing  mechanical  power  tended  to  favor  increased  volume 
amplitudes  at  higher  frequencies  (Figure  1). 

Conclusion:  These  results  in  a  computational  lung  model 
indicate  that  frequency  content  in  MFOV  waveforms  may 
be  tuned  to  minimize  risk  for  VILI.  Our  unique  modeling 
and  optimization  approaches  allow  for  the  selection  of 
patient-specific  MFOV  waveforms,  especially  when 
combined  with  experimental  evidence  to  justify 
physiologically-relevant  objective  emphasis  on  parenchymal 
injury. 
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Figure  1.  Computational  model  of  porcine  lung  (inset). 
Optimal  volume  amplitudes  at  each  frequency  are  shown, 
representing  the  mean  and  standard  deviation  of  the 
waveforms  producing  the  lowest  percentile  of  mechanical 
power  in  healthy  (blue)  and  injured  (red)  couditions. 
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ABSTRACT 


48  High  frequency  oscillatory  ventilation  (HFOV)  relies  on  low  tidal  volumes  cycled  at 

49  supraphysiologic  rates,  producing  fundamentally  different  mechanisms  for  gas  transport  and 

50  exchange  compared  with  conventional  mechanical  ventilation.  Despite  the  appeal  of  using  low 

51  tidal  volumes  to  mitigate  the  risks  of  ventilator-induced  lung  injury  (VILI),  HFOV  has  not 

52  improved  mortality  for  most  clinical  indications.  This  may  be  due  to  non-uniform  and 

53  frequency -dependent  distribution  of  flow  throughout  the  lung.  The  goal  of  this  study  was  to 

54  compare  parenchymal  strain  heterogeneity  during  eucapnic  HFOV  when  using  oscillatory 

55  waveforms  that  consisted  of  either  a  single  discrete  frequency  or  two  simultaneous  frequencies. 

56  We  utilized  on  a  three-dimensional,  anatomically-structured  canine  lung  model  for  simulating 

57  frequency -dependent  ventilation  distribution.  Gas  transport  was  simulated  via  direct  alveolar 

58  ventilation,  advective  mixing  at  bifurcations,  turbulent  and  oscillatory  dispersion,  and  molecular 

59  diffusion.  Volume  amplitudes  at  each  oscillatory  frequency  were  iteratively  optimized  to  attain 

60  eucapnia.  Ventilation  using  single-frequency  HFOV  demonstrated  increasing  heterogeneity  of 

61  acinar  flow  and  CO2  elimination  with  frequency,  for  frequencies  greater  than  the  resonant 

62  frequency.  For  certain  pairs  of  frequencies,  a  linear  combination  of  the  two  corresponding 

63  ventilation  distributions  yielded  reduced  acinar  strain  heterogeneity  compared  to  either  frequency 

64  alone.  Our  model  demonstrates  that  superposition  of  two  simultaneous  oscillatory  frequencies 

65  can  achieve  more  uniform  ventilation  distribution,  and  therefore  lessen  the  potential  for 

66  ventilator-induced  lung  injury,  compared  to  traditional  single-frequency  HFOV. 

67  (219  /  250  words) 

68 


69 


70  NEW  &  NOTEWORTHY 


71  In  this  study  we  simulated  oscillatory  ventilation  with  multiple  simultaneous  frequencies  using  a 

72  computational  lung  model  that  includes  distributed  flow  and  gas  transport.  A  mechanism  of 

73  benefit  was  identified  by  which  ventilation  with  two  simultaneous  frequencies  results  in  reduced 

74  acinar  strain  heterogeneity  compared  to  either  frequency  alone.  This  finding  suggests  the 

75  possibility  of  tuning  the  spectral  content  of  ventilator  waveforms  according  to  patient-specific 

76  mechanical  heterogeneity. 

77  (65/ 75  words) 

78 

79  Keywords:  Computational  models,  Biological  models,  oscillatory  ventilation.  Ventilation 

80  distribution.  Disease  models.  Animal,  Lung  pathology.  Respiratory  mechanics 
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83 

NOMENCLATURE 

84 

Area  of  cross-section  S 

85 

ARDS 

Acute  respiratory  distress  syndrome 

86 

CO2 

Carbon  dioxide 

87 

CT 

Computed  tomography 

88 

d 

Diameter  of  airway  segment 

89 

Deff 

Effective  diffusivity 

90 

^mol 

Molecular  diffusivity 

91 

nlam 

^dis 

Laminar  oscillatory  dispersion  coefficient 

92 

nturb 

^dis 

Turbulent  oscillatory  dispersion  coefficient 

93 

f 

Oscillation  frequency  in  Hz 

94 

fres 

Resonant  frequency  in  Hz 

95 

F 

Fluid  velocity  amplitude  distribution 

96 

G 

Concentration  amplitude  distribution 

97 

HFOV 

High-frequency  oscillatory  ventilation 

98 

i 

Unit  imaginary  number,  equal  to  ^J—l 

99 

In 

Peak  volumetric  strain  in  the  n'^  acinus 

100 

k 

Index  through  oscillation  frequencies 

101 

K 

Number  of  discrete  oscillation  frequencies  in  multi-frequency  waveform 

102 

1 

Length  of  airway  segment 

103 

MFOV 

Multi-frequency  oscillatory  ventilation 

104 

MnCOz 

Molar  flux  of  CO2  for  n*  acinus 

105 

MC02 

Vector  of  acinar  CO2  molar  fluxes 

106 

n 

Index  through  individual  acini 

107 

N 

Total  number  of  acini 

108 

P 

^  atm 

Atmospheric  pressure 

109 

APCO2 

Vector  of  CO2  partial  pressure  differentials  between  acini  and  airway  opening 

no 

r 

Correlation  coefficient 

111 

Jl 

Universal  gas  constant 

112 

M-} 

Real  part  of  the  complexed- valued  argument 

113 

Re 

Reynolds  number 

114 

Recrit 

Critical  Reynolds  number,  transition  between  laminar  and  turbulent  flow 

115 

ndiff 

iv'Y 

Diffusive  transport  resistance 

116 

nmix 

i\'Y 

Diffusive  transport  resistance 

117 

ptot 

iv'Y 

Total  transport  resistance 

118 

Rf 

Transport  resistance  matrix 

119 

t 

Time 

120 

T 

Temperature 

121 

Urms 

Root-mean-square  mean-axial  velocity 

122 

Vd 

Total  dead  space  volume 

123 

I^D,n| 

Magnitude  of  oscillatory  dead  space  volume  delivered  to  n*  acinus 

124 

I^F,n| 

Magnitude  of  oscillatory  fresh  gas  volume  delivered  to  n*  acinus 

125 

^mix 

Oscillatory  mixing  volume 

126 

Vn 

Time- varying  volume  of  the  n*  acinus 

127 

^eg 

Volume  of  airway  segment 

128 

Vj 

Tidal  volume 

129 

Eucapnic  tidal  volume 

130 

l^egl 

Magnitude  of  oscillatory  flow  through  airway  segment 

131 

i4o 

Flow  at  the  airway  opening 

132 

Eucapnic  total  acinar  ventilation 

133 

CO2  elimination  at  the  airway  opening  via  diffusion,  dispersion,  and  mixing 

134 

14CO2 

CO2  elimination  at  acinus  n 

135 

V^^^C02 

CO2  elimination  at  acinus  n  via  diffusion,  dispersion,  and  mixing 

136 

14DAVC02 

CO2  elimination  at  acinus  n  via  direct  acinar  ventilation 

137 

l^totCO^ 

Total  CO2  elimination  at  the  airway  opening 

138 

l^tofCO^ 

Eucapnic  total  CO2  elimination  at  the  airway  opening 

139 

Flow  distribution  similarity  index 

140 

5 

Cross-section  of  airway 

141 

VILI 

Ventilator- induced  lung  injury 

142 

w 

Fluid  velocity 

143 

x,y,z 

Spatial  coordinate  axes 

144 

863 

Conversion  factor  between  partial  pressure  in  mmHg  and  gas  volume  fraction. 

145 

accounting  for  the  difference  between  standard  temperature  and  pressure  dry  air, 

146 

and  body  temperature  and  pressure  air  saturated  with  water  vapor. 

147 

P 

Relative  amplitude  scaling  factor  for  multi-frequency  oscillatory  ventilation 

148 

waveform 

149 

Y 

Concentration  gradient 

150 

^V,n 

Volumetric  strain  at  acinus  n 

151 

0 

Concentration 

152 


Radians  in  semi-circle 


TT 

153  V 

154  cj) 

155  4^ 

156 

157  o) 

158 

159 

160 


Kinematic  viscosity 

Phase  offset  of  sinusoidal  flow  oscillation 

Eucapnic  flow  amplitude  scaling  factor  of  multi-frequency  oscillatory  ventilation 

waveform 

Angular  frequency 


161  INTRODUCTION 

162  High-frequency  oscillatory  ventilation  (HFOV)  relies  on  low  tidal  volumes  cycled  at 

163  supraphysiologic  rates,  which  produces  fundamentally  different  mechanisms  for  gas  transport 

164  compared  with  conventional  mechanical  ventilation.  Despite  the  appeal  of  using  low  tidal 

165  volumes  to  mitigate  the  risks  of  ventilator-induced  lung  injury  (VILI)  in  conditions  such  as  the 

166  acute  respiratory  distress  syndrome  (ARDS),  HFOV  has  not  improved  mortality  in  most  clinical 

167  applications  (8,  22,  23,  38).  The  nonsuperiority  of  HFOV  over  conventional  mechanical 

168  ventilation  (CMV)  may  be  due  to  its  tendency  to  distribute  flow  throughout  the  lung  in  a  non- 

169  uniform  and  frequency -dependent  manner,  especially  in  the  presence  of  lung  tissue  mechanical 

170  heterogeneity  (4,  6).  Reductions  in  prescribed  tidal  volume  during  HFOV  may  therefore  be 

171  counterbalanced  by  increased  strain  heterogeneity  (13).  That  is,  oscillatory  ventilation  at  a  single 

172  discrete  frequency  may  result  in  some  lung  regions  being  underventilated  and  predisposed  to 

173  atelectrauma,  while  other  regions  are  overventilated  and  subject  to  volutrauma  (3).  VILI  has  also 

174  been  associated  with  increased  strain  rate  (26)  and  mechanical  power  (11),  which  depend  on  both 

175  frequency  and  amplitude  of  oscillatory  stretch  and  therefore  may  be  regionally  amplified  during 

176  heterogeneous  ventilation.  Thus,  maldistribution  of  oscillatory  flow  may  increase  the  risk  for 

177  worsening  VILI.  Furthermore  heterogeneous  distribution  of  flow  may  contribute  to  ventilation- 

178  perfusion  mismatching  and  inefficient  gas  exchange  (13). 

179  In  a  previous  study  of  HFOV  in  preterm  lambs  (16),  we  demonstrated  that  oscillation 

1 80  with  multiple  simultaneous  frequencies  resulted  in  improved  gas  exchange  and  lung  recruitment 

181  at  lower  distending  pressures,  compared  to  more  traditional  ‘single-frequency’  oscillatory 

182  ventilation.  We  hypothesized  that  these  short-term  physiologic  improvements  resulted  from  a 

183  more  even  distribution  of  ventilation  throughout  the  heterogeneous  surfactant-deficient  lung.  We 
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184  further  conjectured  that  during  ‘multi- frequency’  oscillatory  ventilation  (MFOV),  flow 

185  oscillations  are  distributed  throughout  the  injured  lung  in  accordance  with  local  mechanical 

186  properties  of  the  parenchyma  and  air- liquid  interface.  Thus  particular  regions  of  the  lung  may  be 

187  preferentially  ventilated  by  some  frequencies,  yet  simultaneously  reject  other  frequencies  (4,  13). 

188  Such  a  mechanical  filtering  effect  during  MFOV  waveform  may  maximize  the  participation  of 

189  various  lung  regions  in  gas  exchange,  despite  different  ‘preferred’  frequencies  within  each 

190  region.  As  a  result,  MFOV  may  produce  more  efficient  CO2  elimination  compared  to  traditional 

191  single-frequency  oscillatory  ventilation,  along  with  a  possible  reduction  in  net  parenchymal 

192  strain  heterogeneity  and  less  potential  for  VILI. 

193  In  this  study,  we  hypothesized  that  oscillatory  ventilation  with  two  simultaneous 

194  frequencies  results  in  more  uniform  parenchymal  strain  distribution  compared  to  oscillation  with 

195  either  frequency  alone,  with  a  reduction  in  peak  strains  required  to  maintain  eucapnic  CO2 

196  elimination.  To  test  this  hypothesis  we  relied  on  our  previously  published  computational  model 

197  of  the  canine  lung  (13),  which  allowed  for  simulated  gas  transport  via  direct  alveolar  ventilation, 

198  mixing  at  bifurcations,  turbulent  and  oscillatory  dispersion,  as  well  as  molecular  diffusion.  We 

199  compared  the  distributions  of  parenchymal  strain  during  traditional  single  frequency  HFOV,  and 

200  then  selected  candidate  pairs  of  frequencies  for  dual-frequency  oscillatory  ventilation  based  on  a 

201  ventilation  distribution  similarity  index.  Finally,  we  assessed  whether  such  a  similarity  index 

202  can  be  an  appropriate  criterion  for  selecting  complementary  frequencies  during  dual-frequency 

203  oscillatory  ventilation. 

204 
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205  METHODS 

206  Model  Overview 

207  A  complete  description  of  our  computational  lung  model  can  be  found  in  Herrmann  et  al. 

208  (13).  Briefly,  the  model  structure  is  based  on  a  three-dimensional  central  airway  tree,  segmented 

209  from  a  computed  tomographic  image  of  a  supine  dog.  Peripheral  airways  with  diameters  smaller 

210  than  2  mm  were  generated  using  a  space-filling  algorithm  (34).  The  model  consisted  of  40,085 

211  cylindrical  airway  segments  in  total,  with  19,932  viscoelastic  acini.  Each  airway  segment  was 

212  characterized  by  a  unique  length  I  and  diameter  d,  with  luminal  volume  given  by: 

(1) 

213  Regional  pleural  pressure  was  assumed  to  vary  with  gravity  according  to  the  hydrostatic  weight 

214  of  the  lung  in  a  supine  orientation  by  0.25  cmH20  per  cm  lung  height  (19,  21),  with  an  average 

215  value  of  -5  cmH20.  Acinar  viscoelasticity  was  determined  according  to  a  second-order 

216  polynomial  function  of  local  transpulmonary  pressure  (4,  6,  18).  Longitudinal  and  shunt 

217  impedances  for  each  airway  segment  accounted  for  pressure  losses  associated  with  oscillatory 

218  flow,  isothermal  gas  compression,  as  well  as  soft  tissue  and  cartilaginous  viscoelasticity  (17). 

219  Flow  oscillations  delivered  at  the  trachea  were  distributed  throughout  the  lung  periphery  using  a 

220  recursive  flow-dividing  scheme,  according  to  the  local  input  impedance  at  each  node  in  the  tree 

221  (6).  Flow  in  each  segment  was  characterized  as  laminar  or  turbulent  according  to  the  Reynolds 

222  number  Re  (3 1 ) : 

Re  =  ^^,  (2) 
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232 
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234 

235 

236 

237 

238 

239 

240 


where  1/rms  is  the  root-mean-square  of  the  time-varying  mean  axial  velocity  and  v  is  kinematic 
viscosity  of  the  gas.  The  transition  from  laminar  to  turbulent  oscillatory  flow  was  assumed  to 
occur  at  the  critical  Reynolds  number  Re^rit  =  30  (14,  31). 

Gas  Transport  Mechanisms 

The  transport  of  CO2  throughout  the  airway  network  was  solved  as  a  steady-state 
problem,  being  time-averaged  throughout  the  oscillatory  ventilation  cycle  and  with  total  CO2 
elimination  equal  to  metabolic  CO2  production  (13).  The  molar  flux  of  CO2  through  each  airway 
segment  was  driven  by  a  concentration  gradient,  and  characterized  by  its  effective  diffusivity 


D. 


eff- 


D 


eff 


K 

'mol  +  y  Oirc/O  . 

,  Re  <  Re^rit 

k=l 

1  -1- 

mol  '  ^dis  ' 

.  Re  ^ 

(3) 


During  laminar  oscillatory  flow  at  multiple  frequencies  (i.e..  Re  <  Recrit)>  l^eff  is  equal  to  the 
sum  of  its  molecular  diffusivity  and  the  laminar  dispersion  coefficients  computed  at 
each  discrete  frequency  (Appendix  A-1).  However  during  turbulent  oscillatory  flow,  Dgff  is 
equal  to  the  sum  of  and  the  frequency-independent  dispersion  coefficient  (9,  35)  of  the 
form: 


=  0.71/RMsd  ■  (4) 

Equation  (4)  is  consistent  with  previous  computational  studies  (9,  30)  based  on  theoretical  (35) 
and  experimental  results  (29).  Oscillatory  flows  were  assumed  to  be  either  fully  laminar  or  fully 
turbulent,  ignoring  transitional  flows.  The  resulting  CO2  transport  resistance  due  to  diffusion 
and  dispersion  is  given  by: 


4JITI 

Tcd^D^ff 


(5) 
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241  where  “R  is  the  universal  gas  constant  and  T  is  absolute  temperature  of  the  gas. 

242  In  our  previous  model  of  convective  gas  transport  (13),  we  described  the  linear  transport 

243  resistance  due  to  gas  mixing  during  repeated  division  and  recombination  of  oscillatory 

244  flow  past  airway  bifurcations  during  oscillation  at  a  single  discrete  frequency.  With  this 

245  approach,  is  characterized  by  the  oscillation  frequency,  /,  and  the  mixing  volume, 

246  which  quantifies  the  volume  of  gas  that  passes  completely  through  an  airway  during  each  period 

247  of  oscillation.  However  during  oscillation  with  multiple  simultaneous  frequencies,  is 

248  computed  over  one  period  of  the  fundamental  frequency  (/o)  of  the  broadband  flow  waveform: 

1/ fo 

j  Ki.it)  dt,  (6) 

0 

249  where  Ki.  is  the  flow  contributing  to  given  by: 

V  .  (t)=  .7^ 

I  0  ,  otherwise 

250  V(.t)  is  the  time-varying  flow  in  an  airway  segment,  l^gg  is  the  lumenal  volume  of  the  airway 

25 1  segment,  and  V (t)  is  the  time  integral  of  flow: 

t 

V(t)  =  j  V(t)  dr  .  (8) 

0 

252  In  practice,  may  be  computed  using  numerical  integration  techniques  (15).  An  example 

253  flow  waveform  through  a  hypothetical  airway  segment  of  volume  70  mL,  along  with  the 

254  corresponding  Kix  ^nd  V^i-^  for  the  segment,  is  shown  in  Figure  1.  For  the  special  case  of 

255  single-frequency  oscillation,  this  numerical  computation  of  Fj^ix  is  consistent  with  the  analytical 

256  expression  for  V^i-^  in  our  previous  model  (13).  The  resulting  gas  transport  resistance  due  to 

257  mixing,  is  then  given  by: 
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The  total  transport  resistance,  in  each  airway  segment  is  given  by  the  parallel  combination 
of  diffusive/dispersive  and  mixing  resistances: 


ntot  _ 
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ndiff  '  nmix 


+ 


(10) 


describes  the  relationship  between  CO2  molar  flux  through  an  airway  segment  and  CO2 
partial  pressure  difference  across  the  airway  segment.  Combining  the  transport  resistances  of  all 
airways  in  the  model  yields  a  system  of  equations  describing  the  relationship  between  acinar  CO2 
partial  pressure  PnC02  and  acinar  CO2  molar  flux  due  to  diffusion,  dispersion,  and  mixing 


M^CO^: 


RjMC02  =  ^PC02  ,  (11) 

where  Rj  is  a  full  matrix  of  transport  resistances,  MCO2  is  a  vector  of  acinar  CO2  molar  fluxes, 
and  APCO2  is  a  vector  of  CO2  partial  pressure  differentials  between  acini  and  airway  opening 
(13).  After  solving  Equation  (11)  to  obtain  MCO2,  acinar  CO2  elimination  due  to  diffusion, 
dispersion,  and  mixing,  V^^^C02,  is  computed  as: 

JiT 

V^^^C02=-^MnC02.  (12) 

^atm 

During  ventilation  with  large  tidal  volumes,  gas  transport  occurs  primarily  via  direct  ventilation, 
whereby  fresh  gas  is  advected  directly  between  the  airway  opening  and  acini.  Under  such 
circumstances,  CO2  elimination  from  the  n*  acinus  due  to  convective  ventilation,  14°'^^C02,  is 
determined  by  the  rate  of  fresh  gas  delivery: 
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273  where  Pn^^2  1®  the  CO2  partial  pressure  of  the  n*  acinus,  863  is  a  conversion  factor  between 

274  partial  pressure  in  mmHg  and  gas  volume  fraction  (13),  and  V-^^nifk)  is  the  fresh  gas  delivered  to 

275  the  n*  acinus  per  oscillation  at  the  /c*  frequency: 

I Vp.n (//c)  I  =  max  (^0,  -  I Vd*  „  | ^  (14) 

276  where  |l4(//c)|  is  the  magnitude  of  oscillatory  flow  delivered  at  the  frequency.  Total  CO2 

277  elimination  for  each  acinus,  l^COz^  is  computed  as  the  sum  of  V-^^^C02  and  V^^^C02.  Total 

278  CO2  elimination  for  the  entire  lung  model,  VtotCOz^  is  then  the  sum  of  I4CO2  across  all  acini: 

N 

l4otC02  =  +  v^^^C02) .  (15) 

n=l 

279  Figure  2  illustrates  which  of  the  modeled  gas  transport  mechanisms  dominates  overall  gas 

280  transport  during  eucapnic  single-frequency  ventilation  between  0.3  Hz  and  26  Hz.  Direct  acinar 

281  ventilation  accounts  for  the  majority  of  total  respiratory  gas  exchange  at  low  frequencies  (i.e., 

282  below  2  Hz),  yet  is  entirely  absent  at  high  frequencies  (i.e.,  above  5  Hz).  The  greatest 

283  contribution  to  the  overall  rate  of  gas  transport  in  individual  airways  derives  from  the  transport 

284  mechanism  describing  bifurcation  mixing  due  to  the  repetitive  division  and  recombination  of 

285  gases.  The  contribution  of  turbulent  mixing  to  gas  transport  in  individual  airways  increases  at 

286  higher  frequencies,  dominating  in  the  larger  central  airways  during  delivery  of  small  volume 

287  amplitudes  at  high  instantaneous  flow  rates.  This  behavior  is  consistent  with  the  analysis  of 

288  bifurcating  flow  regimes  by  Jan  et  al.  (14)  as  well  as  the  numerical  simulations  performed  by 

289  Choi  et  al.  (5). 
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290 


Comparison  of  Single -Frequency  and  Dual-Frequency  Oscillation 


291  To  identify  plausible  candidate  frequency  pairs  for  dual-frequency  simulations,  each 

292  combination  of  two  distinct  oscillatory  frequencies  was  considered  by  taking  the  inner  product 

293  between  their  respective  normalized  acinar  flow  distributions: 

N 

n=l 

Sv(k.f2)  = 

Xlv-nC/or  ^|r„(/2)f 

^n=l  ^n=l 

294  where  — the  value  of  the  inner  product — quantifies  the  ventilation  distribution  similarity. 

295  Thus  Sy  =  1  implies  two  identical  distributions,  whereas  =  0  implies  two  orthogonal 

296  distributions.  Figure  3  shows  a  schematic  lung  with  three  mechanically  heterogeneous 

297  compartments  as  a  simplified  example  of  quantifying  ventilation  distribution  similarity. 

298  Figure  4A  shows  the  contours  of  for  each  permutation  of  two  frequencies  from  0.2 

299  Hz  to  40  Hz,  spaced  over  0.2  Hz  intervals.  Identifying  pairs  of  frequencies  with  the  lowest 

300  degree  of  similarity  (i.e.,  the  smallest  values  of  Sy)  reveals  the  most  complementary  pairs,  i.e., 

301  those  pairs  which  are  most  able  to  compensate  for  each  other’s  respective  underventilated  and 

302  overventilated  regions.  Alternatively,  the  Pearson  correlation  coefficient  (r)  between  normalized 

303  acinar  flow  distributions  for  each  frequency  pair  can  also  be  used  to  identify  complementary 

304  distributions.  Positive  values  of  r  correspond  to  similar  distributions,  such  that  acini  tend  to 

305  receive  either  above-average  or  below-average  amounts  of  flow  at  both  frequencies.  By  contrast, 

306  negative  values  of  r  correspond  to  dissimilar  distributions,  such  that  acini  tend  to  receive  above- 

307  average  flow  at  one  frequency  yet  below-average  at  the  other.  Figure  4B  shows  the  contours  of 

308  r  for  each  permutation  of  two  frequencies  between  0.2  Hz  and  40  Hz.  Identifying  pairs  of 
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frequencies  with  the  most  negative  correlation  coefficients  reveals  the  most  complementary 
pairs. 

Based  on  the  contour  maps  shown  in  Figure  4,  we  selected  a  pair  of  frequencies 
representing  a  low  value  of  (i.e.,  dissimilar  ventilation  distributions)  for  simulation  of  gas 
transport.  We  selected  the  local  minimum  of  the  contour  map  corresponding  to  12  Hz  and  26 
Hz.  For  comparison,  another  pair  of  frequencies  with  a  high  value  of  was  also  selected  for 
simulation,  at  0.3  Hz  and  12  Hz.  For  each  pair  of  frequencies,  multiple  waveforms  were 
constructed  to  define  flow  at  the  airway  opening  (l4o)  according  to  the  following  template: 

14o(t)  =  ^euclP  sin(27r/it)  +  (1  -  p)  sin(27r/2t  +  0)] ,  (17) 


where  P  is  a.  scaling  factor  adjusting  the  relative  contributions  of  /i  versus  /2  such  that  0  <  P  < 
1,  0  is  the  relative  phase  offset  between  the  flow  oscillations  such  that  0  <  0  <  2n,  and 
^euc  >  0  is  a  scaling  factor  used  to  adjust  the  flow  amplitudes  delivered  by  the  MFOV 
waveform.  For  any  combination  of  /i,  /2,  P,  and  0,  the  corresponding  value  of  H^guc  was 
iteratively  optimized  to  attain  eucapnic  CO2  elimination,  which  was  approximated  according  to 
the  assumed  weight  of  the  animal  (13).  For  simplicity,  the  value  of  0  was  constrained  to  be  zero 
for  all  simulations. 

Following  the  determination  of  eucapnic  scaling,  the  distribution  of  peak  acinar 
volumetric  strain  (sy  n)  was  computed  by  calculating  the  range  of  time- varying  acinar  volume, 
14  (t),  divided  by  the  minimum  value: 


max[F^(t)]  -  min[F^(t)] 
min[F„(t)] 


(18) 
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327  The  distribution  of  acinar  root-mean- square  volume  (Vrms.ti)  was  also  used  to  assess  distributed 

328  ventilation  heterogeneity: 

I  i7^ 

W,n=  /o  I  [Vn(t)-V,^Vdt  .  (19) 

^|  0 

329  where  is  time-averaged  acinar  volume  given  by: 

'^1  fo 

K=fol  Vnmt.  (20) 

0 

330  The  values  of  fRMS,n  ^^d  in  Equations  (19)  and  (20)  were  numerically  approximated  using 

331  trapezoidal  integration  (15). 

332  Simulations 

333  The  recursive  algorithms  for  computing  and  storing  impedances  and  determining 

334  advective  flow  division  throughout  the  tree  were  written  in  C-t-i-.  Computation  time  for 

335  generating  advective  flow  distribution  at  each  distinct  /  was  ~3  s  on  a  64-bit  computer  with  an 

336  Intel  Core  17-950  processor  operating  at  3.07  GHz  with  12  GB  RAM.  The  gas  transport  model 

337  was  written  and  executed  in  MATLAB  (version  R2016a,  The  Mathworks,  Natick,  MA).  Total 

338  computation  time  for  each  distinct  dual-frequency  simulation  was  ~70  s.^ 

339 

340 


*  A  sample  MATLAB  script  for  computing  the  gas  transport  model  is  available  from  the  authors 
upon  request. 
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341 


RESULTS 


342  Two  pairs  of  frequencies  were  selected  for  gas  transport  simulations  based  on  ventilation 

343  distribution  similarity  and  correlation  coefficient  (r),  as  indicated  by  the  labeled  points  on 

344  the  contour  maps  of  Figure  4.  The  frequency  pair  with  similar  distributions  was  chosen  at  0.3  Hz 

345  and  12  Hz,  with  =  0.98  and  r  =  +0.69.  The  pair  with  dissimilar  distributions  was  chosen  at 

346  12  Hz  and  26  Hz,  with  =  0.85  and  r  =  —0.64.  Figure  5  shows  the  anatomic  distributions  of 

347  acini  receiving  relatively  more  or  less  oscillatory  flow  at  0.3  Hz,  12  Hz,  and  26  Hz  compared  to  a 

348  perfectly  uniform  distribution.  Acinar  flow  magnitudes  were  normalized  relative  to  their 

349  corresponding  values  for  a  symmetric,  homogeneous  lung  with  rigid  airway  walls  and  no  gas 

350  compression  (6).  Figure  6  shows  scatter  plots  comparing  normalized  ventilation  distributions 

351  between  frequencies  for  both  selected  frequency  pairs.  The  dissimilar  distributions  at  12  Hz  and 

352  26  Hz  exhibit  a  majority  of  acini  (72.7  %)  that  are  preferentially  ventilated  by  only  one 

353  frequency  or  the  other,  whereas  the  similar  distributions  at  0.3  Hz  and  12  Hz  exhibit  only  a 

354  minority  of  acini  (17.3  %)  with  the  same  complementary  behavior. 

355  Figure  7  shows  the  required  eucapnic  flow  scaling  for  dual-frequency  waveforms 

356  constructed  from  various  linear  combinations  of  each  frequency  in  the  selected  frequency  pairs, 

357  according  to  Equation  (17).  For  all  values  of  p,  values  of  T'euc  were  iteratively  adjusted  to 

358  achieve  total  CO2  elimination  within  0.3  %  error  of  the  predicted  metabolic  CO2  production  of 

359  1.93  •  10“^  L  for  a  25  kg  dog.  Figure  8  shows  the  resulting  distributions  of  acinar  peak 

360  volumetric  strain  (£v,ti)>  Figure  9  shows  distributions  of  acinar  root- mean- square  volume 

361  (Irms.ti)-  Comparing  distributed  ventilation  using  oscillations  at  only  a  single  frequency  (i.e.,  at 

362  =  0  or  =  1  only),  0.3  Hz  ventilation  produced  the  least  heterogeneity  in  both  strain  and 

363  Frmsj  but  also  the  largest  magnitudes  of  strain,  compared  to  either  12  Hz  or  26  Hz.  For  example. 
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364  the  range  of  reported  as  median  (minimum,  maximum),  was  15.5  %  (12.3  %,  19.4  %) 

365  during  single-frequency  ventilation  at  0.3  Hz,  compared  to  2.8  %  (6.5  %,  2.1  %)  at  12  Hz  or  1.6 

366  %  (0.5  %,  3.4  %)  at  26  Hz.  For  the  dual-frequency  waveforms  constructed  using  0.3  Hz  and  12 

367  Hz,  no  improvement  in  ventilation  heterogeneity  was  observed  for  any  linear  combination  of  the 

368  two  frequencies,  compared  to  either  single  frequency  alone  (i.e.,  0  <  p  <  1  compared  to  p  =  0 

369  or  p  —  1).  However  for  waveforms  using  12  Hz  and  26  Hz  oscillations,  an  optimal  value  of 

370  p  =  0.4  was  found  to  provide  the  greatest  reduction  in  the  difference  between  the  maximum  and 

371  minimum  values  of  normalized  acinar  strain  and  Vrms-  Specifically  at  p  —  0.4,  the  range  of 

372  normalized  strain  was  reduced  by  50.8  %  compared  to  26  Hz  (P  —  0)  and  37.2  %  compared  to 

373  12Hz(;g  =  l). 

374 
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375  DISCUSSION 

376  In  this  study,  we  extended  an  existing  computational  model  of  distributed  flow  and  gas 

377  transport  (13)  to  simulate  ventilation  heterogeneity  during  oscillatory  ventilation  with  either  a 

378  single  frequency  or  two  frequencies  delivered  simultaneously.  By  choosing  pairs  of  frequencies 

379  with  complementary  ventilation  distributions,  we  demonstrated  that  regional  underventilation  at 

380  one  frequency  may  be  compensated  by  regional  overventilation  at  the  other  frequency.  Thus,  an 

381  ‘optimal’  linear  combination  of  two  frequencies  can  be  identified  that  yields  the  least  amount  of 

382  parenchymal  stain  heterogeneity  compared  to  either  single  frequency  alone  (Figures  8B  and  9B). 

383  However,  for  a  pair  of  oscillatory  frequencies  yielding  similar  ventilation  distributions,  no  linear 

384  combination  was  found  to  reduce  ventilation  heterogeneity  (Figures  8 A  and  9A).  This  finding 

385  suggests  that  it  is  possible  to  optimize  the  spectral  content  of  oscillatory  waveforms  to  achieve 

386  substantial  reductions  in  ventilation  heterogeneity  and  parenchymal  strain. 

387  In  these  simulations  of  eucapnic  oscillatory  ventilation  in  a  canine  lung,  low  frequencies 

388  produced  the  lowest  degree  of  ventilation  heterogeneity,  yet  the  largest  magnitudes  of  regional 

389  strain.  However,  increasing  oscillatory  frequency  resulted  in  reduced  regional  strain  but 

390  increased  strain  heterogeneity,  especially  for  frequencies  greater  than  the  resonant  frequency  of 

391  this  model,  f^-es  —  7.5  Hz  (13).  The  resonant  frequency  indicates  the  transition  from  flows 

392  dominated  by  elastic  forces  (/  <  /res)  to  flows  dominated  by  inertial  forces  (/  >  /res)  (4,  6). 

393  Thus  high-frequency  oscillatory  ventilation  (HFOV)  may  reduce  parenchymal  strain,  yet 

394  increase  both  strain  rate  and  strain  heterogeneity  compared  to  conventional  mechanical 

395  ventilation  (20,  25,  26).  This  trade-off  may  in  part  explain  the  nonsuperiority  of  HFOV  over 

396  conventional  mechanical  ventilation  in  patients  with  mechanically  heterogeneous  lungs  (8,  22, 

397  38). 
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398  Our  computational  model  highlights  a  means  by  which  both  parenchymal  strain  and 

399  strain  heterogeneity  may  be  minimized.  If  multiple  frequencies  are  delivered  simultaneously,  the 

400  respective  contributions  of  each  frequency  to  the  combined  distribution  of  parenchymal  strain 

401  may  be  tuned  towards  a  particular  objective,  e.g.,  minimizing  strain  heterogeneity.  By  extending 

402  this  concept  of  dual-frequency  oscillation  to  the  more  general  case  of  three  of  more  frequencies 

403  (16),  one  may  further  exploit  frequency-dependent  ventilation  distribution  to  mitigate  adverse 

404  consequences  of  parenchymal  strain  heterogeneity.  This  proposed  mechanism  of  benefit  may  in 

405  part  explain  improvements  in  oxygenation  and  gas  exchange  reported  during  experimental  use  of 

406  multi-frequency  oscillatory  ventilation  (MFOV)  in  surfactant-deficient  preterm  lungs  (16). 

407  These  findings  are  also  relevant  to  high-frequency  percussive  ventilation  (HFPV),  which 

408  involves  high  frequency  oscillations  superimposed  on  conventional  mechanical  ventilation 

409  waveforms  (1,  7,  28).  Large  reductions  in  parenchymal  strain  at  lower  frequencies  may  be 

410  achieved  during  HFPV  by  supplementary  gas  transport  using  a  relatively  small  contribution  of 

411  superimposed  high  frequency  oscillation. 

412  Two  criteria  were  proposed  to  evaluate  candidate  frequency  pairs  in  this  study:  the 

413  ventilation  distribution  similarity  (cSj^)  and  the  correlation  coefficient  (r).  As  shown  in  Figure  4, 

414  and  r  demonstrate  close  agreement  with  each  other  over  the  range  of  frequencies  from  0.1  Hz 

415  to  30  Hz,  suggesting  that  the  outcome  of  candidate  frequency  pair  selection  may  not  be  specific 

416  to  the  particular  criterion  used  for  evaluation. 

417  Although  peak  acinar  volumetric  strain  and  specific  ventilation  are  typically  used  to 

418  assess  regional  ventilation  heterogeneity  during  tidal  breathing  (2,  27),  these  metrics  may  not  be 

419  directly  applicable  during  oscillation  with  two  or  more  frequencies,  since  the  corresponding 

420  regional  volume  fluctuations  will  not  oscillate  between  a  consistent  minimum  and  maximum. 
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421  Rather,  parenchymal  strain  may  exhibit  various  local  minima  and  maxima.  Thus  ventilation 

422  heterogeneity  during  multi-frequency  oscillation  may  be  more  appropriately  characterized  by  the 

423  root-mean-square  of  regional  volume  fluctuations  (Rrms)  (16).  Vrms  describes  the  variance  of 

424  regional  distension  about  a  mean  lung  volume,  which  is  proportional  to  volume  amplitude  during 

425  single-frequency  oscillatory  ventilation,  yet  also  accounts  for  oscillation  at  multiple  distinct 

426  frequencies.  In  Figures  8  and  9  the  distributions  of  peak  acinar  volumetric  strain  (^v.n) 

427  Frms.ti  exhibited  close  agreement,  suggesting  that  the  analysis  of  ventilation  heterogeneity  is  not 

428  sensitive  to  the  choice  of  either  metric. 

429  Limitations 

430  The  modeling  assumptions  and  validation  for  the  single-frequency  gas  transport  model 

431  have  been  presented  in  our  previous  study  (13).  Additional  comparison  to  distributed  strain 

432  measured  on  the  pleural  surface  of  excised  canine  lungs  is  provided  in  Appendix  A-2.  Our 

433  currently  proposed  extension  to  accommodate  two  frequencies  of  oscillations  simultaneously  has 

434  been  derived  using  a  similar  conceptual  framework.  However  to  our  knowledge,  no  robust 

435  experimental  data  is  available  to  validate  these  model  simulations.  Specific  or  quantitative 

436  predictions  based  on  this  model  should  therefore  be  treated  with  appropriate  caution. 

437  Particularly,  our  simulations  assume  that  the  mechanical  behavior  of  the  lung  is  linear  and  time- 

438  invariant.  Thus  superposition  of  two  sinusoidal  oscillations  at  the  airway  opening  results  in  a 

439  simple  linear  combination  of  the  respective  flow  distributions  at  each  frequency.  However, 

440  nonlinear  mechanical  properties  of  an  actual  mammalian  lung  may  lead  to  unpredictable 

441  distortion  and  interaction  of  superimposed  oscillations  (32,  33,  39).  Experimental  imaging 

442  techniques  such  as  frequency -selective  computed  tomography  (12)  or  optical  coherence 
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443  tomography  (24)  may  be  used  to  further  assess  the  extent  of  nonlinear  mechanical  behavior 

444  during  multi-frequency  oscillation  in  situ. 

445  This  particular  analysis  of  simulated  multi-frequency  oscillation  examined  only  two 

446  specific  pairs  of  frequencies,  to  establish  the  feasibility  of  using  a  multi-frequency  approach  to 

447  reduce  parenchymal  strain  heterogeneity  compared  to  single-frequency  oscillation  in  a 

448  computational  model.  Furthermore  the  simulations  were  performed  using  a  model  of  healthy 

449  canine  lungs,  with  all  mechanical  heterogeneity  derived  only  from  gravitational  gradients  in 

450  transpulmonary  pressure  and  airway  network  asymmetries.  Nonetheless  this  healthy  model 

451  exhibited  substantial  ventilation  heterogeneity  at  high  frequencies  (10,  20),  indicating  that 

452  exploiting  frequency-dependent  ventilation  heterogeneity  may  be  possible  even  within  lungs 

453  with  homogeneous  tissue  compliance.  The  particular  distribution  of  regional  strain  may  be 

454  appreciably  altered  under  injured  conditions  (13),  yet  the  rationale  for  applying  a  multi- 

455  frequency  approach  remains  the  same  (16).  The  results  of  these  simulations  warrant  further 

456  investigation,  involving  optimization  of  the  spectral  content  in  oscillatory  waveforms,  or  using  a 

457  larger  number  of  distinct  frequencies.  Such  an  optimization  procedure  may  involve  an 

458  evaluation  of  the  ventilation  distribution  similarity,  similar  to  that  presented  in  Figure  4. 

459  However  the  prospect  of  translating  these  computational  modeling  results  into  a  practical  in  vivo 

460  approach  for  identification  of  optimal  frequencies  presents  several  considerable  technical 

461  challenges,  not  the  least  of  which  is  non-invasive  quantification  of  frequency-dependent 

462  ventilation  heterogeneity  (10,  27). 

463 
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464 


CONCLUSION 


465  Although  high  frequency  oscillatory  ventilation  achieves  CO2  elimination  with  reduced 

466  parenchymal  strain  compared  to  conventional  mechanical  ventilation,  increased  strain 

467  heterogeneity  during  HFOV  may  be  injurious  to  the  lung.  Heterogeneity  of  regional  mechanical 

468  properties  and  airway  network  asymmetry  result  in  frequency-dependent  distributions  of 

469  oscillatory  flow,  especially  at  high  frequencies.  Such  frequency-dependent  heterogeneity  in 

470  ventilation  may  be  exploited  by  combining  multiple  simultaneous  frequencies  of  oscillation. 

471  Dual-frequency  oscillatory  ventilation  may  provide  more  uniform  ventilation  throughout  the 

472  heterogeneous  lung,  such  that  regional  maldistribution  of  parenchymal  strain  can  be  minimized 

473  by  using  oscillatory  frequencies  with  dissimilar  ventilation  distributions.  This  unique  oscillatory 

474  modality  may  thus  have  utility  for  the  treatment  of  patients  with  heterogeneous  lung  injury,  by 

475  minimizing  the  extremes  of  cyclic  parenchymal  strain  and  reducing  the  risk  of  ventilator- induced 

476  lung  injury. 
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507  FIGURE  LEGENDS 

508  Figure  1.  Illustration  of  method  for  computing  Fj^ix  during  ventilation  with  multiple  frequencies 

509  of  oscillation.  In  this  example,  two  frequencies  are  used — 3  Hz  and  7  Hz — such  that  /o  =  1  Hz. 

510  The  shading  indicates  where  the  conditions  of  F(t)  >0  (light  and  dark)  and  F(t)  >  l^gg  (dark) 

511  are  satisfied.  Fj^ix  can  be  calculated  from  the  sum  of  flow  integrals  indicated  by  dark  shading  for 

512  which  both  conditions  are  satisfied.  Equivalently,  Fj^ix  can  be  calculated  from  the  sum  of  the 

513  lengths  of  the  dark  lines  denoting  the  local  maxima  of  F (t)  —  l^eg- 

514  Figure  2.  Airway  segments  and  acini  labeled  according  to  gas  transport  mechanisms.  (A) 

515  Decision  trees  used  for  labeling  airway  segments  according  to  the  dominant  gas  transport 

516  mechanism  contributing  to  the  overall  rate  of  gas  transport  within  each  individual  airway.  Acini 

517  are  labeled  if  fresh  gas  is  advected  via  bulk  flow  directly  into  the  acinar  space.  (B)  Airways  and 

518  acini  colored  according  to  gas  transport  mechanisms  during  eucapnic  ventilation  at  various 

519  frequencies  using  the  tidal  volumes  indicated  in  mL.  Acini  that  do  not  receive  fresh  gas 

520  ventilation  are  not  shown.  Gravitational  direction  (g)  is  indicated  for  a  supine  orientation. 

521  Figure  3.  Two  examples  of  ventilation  distribution  similarity  (cSj/)  and  correlation  coefficient  (r) 

522  in  a  schematic  lung  with  three  mechanically  heterogeneous  compartments.  The  frequency- 

523  dependent  fraction  of  oscillatory  flow  received  by  each  compartment  is  indicated  for  two 

524  frequencies  and  /2.  In  the  case  of  similar  distributions  (A),  the  central  compartment  receives 

525  the  greatest  proportion  of  oscillatory  flow  at  both  frequencies,  and  the  corresponding  value  of  Sy 

526  is  close  to  unity  whereas  r  is  close  to  +1.  In  the  case  of  dissimilar  distributions  (B),  the 

527  compartments  underventilated  at  are  overventilated  at  /2,  and  vice  versa — thus  the  value  of 

528  is  close  to  zero  and  r  is  close  to  —1.  See  text  for  definitions  of  and  r. 
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529  Figure  4.  Contours  of  (A)  ventilation  distribution  similarity  (<5j/),  and  (B)  correlation  coefficient 

530  (r),  in  the  canine  lung  model  for  pairs  of  frequencies  between  0. 1  and  40  Hz.  The  line  of  identity 

531  (dotted  line)  represents  comparisons  between  each  frequency  and  itself,  such  that  =  1 

532  and  =  +1.  The  contours  of  are  symmetric  about  the  line  of  identity,  such  that 

533  =  Sy(f2,fi)  and  r(/i,/2)  =  r(/2,/i).  Labeled  points  correspond  to  frequency  pairs 

534  selected  for  further  examination:  12  Hz  and  26  Hz  (circle),  0.3  Hz  and  12  Hz  (triangle).  See  text 

535  for  definitions  of  and  r. 

536  Figure  5.  Acini  labeled  according  to  relative  share  of  delivered  oscillatory  flow  during  0.3  Hz, 

537  12  Hz,  and  26  Hz  ventilation.  Acinar  flow  is  normalized  by  the  flow  amplitude  delivered  at  the 

538  trachea  divided  by  the  total  number  of  acini,  such  that  values  above  or  below  1  represent  acini 

539  receiving  more  (red)  or  less  (blue)  than  their  respective  shares  of  uniformly  distributed  flow, 

540  respectively.  Note  that  color  bars  are  scaled  differently  for  each  frequency  to  emphasize  the 

541  regional  distributions  of  relative  over-  and  under- ventilation.  Various  view  angles  are  provided, 

542  with  gravitational  direction  (g)  indicated  for  a  supine  orientation  (directed  downwards,  out  of  the 

543  page,  or  upwards). 

544  Figure  6.  Scatter  plots  comparing  normalized  acinar  flow  distributions  at  two  frequencies,  where 

545  each  point  represents  a  single  acinus,  and  horizontal  or  vertical  position  correspond  to 

546  normalized  flow  amplitude  in  each  acinus  at  either  frequency.  Two  examples  are  shown:  (A) 

547  similar  distributions  at  0.3  Hz  vs.  12  Hz  with  =  0.98  and  r  =  -1-0.69  (left),  and  (B)  dissimilar 

548  distributions  at  12  Hz  vs.  26  Hz  with  =  0.85  and  r  =  —0.64  (right).  The  four  quadrants 

549  delineate  acini  that  receive  either  above-average  (i.e.,  greater  than  1.0)  or  below-average  (i.e., 

550  less  than  1.0)  amounts  of  oscillatory  flow  at  one  or  both  frequencies. 
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551  Figure  7.  Eucapnic  flow  scaling  (H^euc)  for  dual-frequency  waveforms  constructed  according  to 

552  Equation  (17)  using  0  <  <  1  to  scale  the  contribution  of  each  frequency,  using  oscillations  at 

553  either  0.3  Hz  and  12  Hz  (A),  or  12  Hz  and  26  Hz  (B).  Each  point  represents  a  single  simulation. 

554  Corresponding  volume  waveforms  at  the  airway  opening  are  shown  for  each  frequency  pair  at 

555  select  values  of  p. 

556  Eigure  8.  Distributions  of  acinar  peak  volumetric  strain  (£v)  for  dual-frequency  waveforms 

557  constructed  according  to  Equation  (17)  using  0  <  P  <  1  to  scale  the  contribution  of  each 

558  frequency,  using  oscillations  at  either  0.3  Hz  and  12  Hz  (A),  or  12  Hz  and  26  Hz  (B). 

559  Distributions  are  represented  by  the  median  (solid  line)  and  full  range  between  minimum  and 

560  maximum  value  (shaded  region).  Normalized  distributions  are  shown  with  the  mean  value 

561  normalized  to  unity  (dotted  line). 

562  Eigure  9.  Distributions  of  acinar  root-mean-square  volume  (Erms)  for  dual-frequency  waveforms 

563  constructed  according  to  Equation  (17)  using  0  <  P  <  1  to  scale  the  contribution  of  each 

564  frequency,  using  oscillations  at  either  0.3  Hz  and  12  Hz  (A),  or  12  Hz  and  26  Hz  (B). 

565  Distributions  are  represented  by  the  median  (solid  line)  and  full  range  between  minimum  and 

566  maximum  value  (shaded  region).  Normalized  distributions  are  shown  with  the  mean  value 

567  normalized  to  unity  (dotted  line). 
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568  Figure  A- 1.  Comparison  between  experimental  and  simulated  acinar  area  strains.  Experimental 

569  data  for  pleural  surface  area  strain  reported  by  Lehr  et  al.  (20)  is  indicated  by  mean  (black 

570  circles)  and  standard  deviation  (error  bars)  at  1  Hz,  15  Hz,  and  30  Hz.  Equivalent  area  strain  for 

571  our  model  assuming  isotropic  deformation  during  simulated  ventilation  is  shown  as  mean 

572  (dashed  grey  line)  and  standard  deviation  (solid  grey  lines).  (A)  Area  strain  amplitude  is 

573  reported  during  oscillation  using  a  constant  ratio  of  7.6  %  oscillatory  volume  amplitude  relative 

574  to  total  lung  volume  (see  text  for  additional  experimental  details).  (B)  Area  strain  phase  is 

575  reported  relative  to  volume  phase  measured  at  the  airway  opening. 

576 

577 
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578  APPENDIX 


579 

580  A-1.  Laminar  Oscillatory  Dispersion 

581  In  this  analysis  of  laminar  oscillatory  dispersion  during  flow  with  multiple  simultaneous 

582  frequencies,  we  follow  the  reasoning  and  nomenclature  of  Watson  (37).  Consider  multi- 

583  frequency  oscillatory  flow  in  a  tube  with  velocity  profile: 

K 

w(x,y,t)  =  ^  iRe{Ffc(x, ,  (A-1) 

k  =  l 

584  where  k  denotes  the  angular  oscillation  frequencies  ro/j  =  2nfi^  from  1  to  K,  j  is  the  unit 

585  imaginary  number,  and  =  Fi^(x,y)  is  velocity  amplitude  distribution  over  the  cross-section  in 

586  the  x-y  plane,  which  is  periodic  in  time  t.  The  operator  denotes  the  real  part  of  the 

587  complex-valued  argument  within  the  brackets.  The  concentration  profile  of  a  contaminant 

588  species  contained  in  the  carrier  gas  will  then  be  of  the  form  (37): 

K 

6{x,y,z,  t)  =  -yz  +  ^  iRe{7Gfc(x,y)e^'"fc'']  ,  (A-2) 

k=l 

589  where  z  is  the  axial  direction,  y  is  a  constant  corresponding  to  a  linear  concentration  gradient  in 

590  the  axial  direction,  and  =  G/j(x,  y)  is  the  modulated  amplitude  distribution  over  the  cross- 

591  section  in  the  x-y  plane  and  is  periodic  in  time  t  with  angular  frequency  such  that  the 

592  concentration  amplitude  distribution  is  given  by  yG^..  The  rate  of  molar  flux  of  the  contaminant 

593  has  advective  and  diffusive  terms  integrated  over  the  cross-section  S  is  given  by  (37): 

rr  r  ^6^ 

JJ  -  Dmoi  dxdy ,  (A-3) 

which  expands  to: 
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594 

595 

596 

597 


598 

599 

600 


601 


602 

603 

604 


yOmol  +  ^■1'^ 

+  l)'|^(Cie«'  +  G^e-«‘)| 

V/c  =  l  J 


(A-4) 


dxdy , 


where  the  asterisks  denote  complex  conjugates.  The  expansion  of  the  product  produces 
oscillatory  terms  (i.e.  terms  preceding  or  proceeding  = 

I  ^  steady  terms  (i.e.  terms  preceding  =  1).  The  mean  rate 

of  flux,  averaged  over  time,  is: 


£  YDmoi  +  Ir^iF^G*^  + 


dxdy , 


(A-5) 


wherein  all  oscillatory  terms  are  reduced  to  zero  in  the  time-averaging,  and  only  the  steady  terms 
remain  (37).  The  effective  diffusivity  is: 

£^eff  =  ^moi(l  +  Fd)  >  (A-6) 

where  is  the  relative  increase  in  diffusivity  compared  to  molecular  diffusivity,  given  by: 


Rn  = 


which  rearranges  to: 


4a 


K 

[f  y + fi'Gj 

molAs  JJs 


dxdy , 


(A-7) 


K 

k.  —  1 


(A-8) 


where  i4s  is  the  area  of  the  cross-section  S,  and  the  term  inside  the  summation  is  identical  to  the 
relative  increase  in  diffusivity  for  any  individual  frequency  of  oscillation  (37).  Thus  the 
effective  diffusivity  is  given  by: 
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(A-9) 


K 

£>eff  =  £>mol  +  ^£>d1s”(//c)- 
k=l 

605  Equation  (A-9)  relies  on  the  reduction  of  time-averaged  oscillatory  terms  to  zero  over  a  long 

606  period  of  oscillation,  i.e.,  the  elimination  of  and  terms  in  the  transition  from 

607  Equation  (A-4)  to  Equation  (A-5).  Therefore,  this  model  of  dispersion  during  simultaneous 

608  multiple  frequencies  may  not  apply  when  the  beat  frequency  (i.e.,  the  difference  between  any 

609  two  frequencies)  is  small. 

610 

611  A-2.  Comparison  with  Experimental  Strain  Measurements 

612  In  our  previous  work  (13),  we  demonstrated  that  the  tidal  volume  required  for  eucapnic 

613  ventilation  during  single-frequency  oscillatory  ventilation  was  consistent  with  the  experimental 

614  findings  of  Venegas  et  al.  (36).  The  distribution  of  volumetric  strain  within  lungs  in  vivo  is 

615  difficult  to  quantify  experimentally,  especially  during  rapid  motion  such  as  during  high- 

616  frequency  oscillatory  ventilation.  Eehr  et  al.  (20)  used  stroboscopic  photography  to  capture 

617  time- varying  motion  of  grid  points  marked  on  the  pleural  surface  of  an  excised  canine  lung. 

618  They  reported  area  strain  measurements  during  oscillations  at  1  Hz,  15  Hz,  and  30  Hz.  We  may 

619  compare  our  simulated  distributions  of  volumetric  strain  to  the  area  strain  distributions  reported 

620  by  Eehr  et  al.  by  first  adjusting  the  simulated  volume  oscillations  to  achieve  the  same  ratio  of 

621  volume  amplitude  to  total  lung  volume.  Eehr  et  al.  reported  using  100  mE  volume  oscillations 

622  applied  to  a  lung  of  total  volume  1320  mE,  expecting  volumetric  strains  of  approximately  7.6  %. 

623  Assuming  isotropic  expansion  and  contraction,  they  expected  the  corresponding  area  strains  to 

624  be  two-thirds  of  the  volumetric  strain,  or  approximately  5.1  %.  Eigure  A-1  shows  the  mean  and 

625  standard  deviation  of  pleural  surface  area  strain  amplitudes  and  phases  reported  by  Eehr  et  ah,  as 
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626  well  as  the  equivalent  acinar  area  strain  amplitudes  and  phases  computed  from  our  simulations 

627  (assuming  isotropic  deformation).  Both  acinar  strain  amplitude  heterogeneity  and  phase 

628  heterogeneity  increase  at  high  oscillation  frequencies  above  resonant  frequency  (7.5  Hz  in  our 

629  model).  This  trend  is  also  consistent  with  increasing  heterogeneity  in  alveolar  pressure 

630  amplitudes  and  phases  measured  by  Fredberg  et  al.  (10),  although  pressure  measurements  alone 

63 1  cannot  be  used  to  infer  volumetric  strain. 

632 
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Herrmann  J,  Tawhai  MH,  Kaczka  DW.  Regional  gas  transport 
in  the  heterogeneous  lung  during  oscillatory  ventilation.  J  Appl 
Physiol  121:  1306-1318,  2016.  First  published  October  7,  2016; 
doi:10.1152/japplphysiol.00097.2016. — Regional  ventilation  in  the 
injured  lung  is  heterogeneous  and  frequency  dependent,  making  it 
difficult  to  predict  how  an  oscillatory  flow  waveform  at  a  specified 
frequency  will  be  distributed  throughout  the  periphery.  To  predict  the 
impact  of  mechanical  heterogeneity  on  regional  ventilation  distribu¬ 
tion  and  gas  transport,  we  developed  a  computational  model  of 
distributed  gas  flow  and  CO2  elimination  during  oscillatory  ventila¬ 
tion  from  0.1  to  30  Hz.  The  model  consists  of  a  three-dimensional 
airway  network  of  a  canine  lung,  with  heterogeneous  parenchymal 
tissues  to  mimic  effects  of  gravity  and  injury.  Model  CO2  elimination 
during  single  frequency  oscillation  was  validated  against  previously 
published  experimental  data  (Venegas  JG,  Hales  CA,  Strieder  DJ,  J 
Appl  Physiol  60:  1025-1030,  1986).  Simulations  of  gas  transport 
demonstrated  a  critical  transition  in  flow  distribution  at  the  resonant 
frequency,  where  the  reactive  components  of  mechanical  impedance 
due  to  airway  inertia  and  parenchymal  elastance  were  equal.  For 
frequencies  above  resonance,  the  distribution  of  ventilation  became 
spatially  clustered  and  frequency  dependent.  These  results  highlight 
the  importance  of  oscillatory  frequency  in  managing  the  regional 
distribution  of  ventilation  and  gas  exchange  in  the  heterogeneous 
lung. 

high-frequency  oscillatory  ventilation;  gas  exchange;  lung  injury; 
respiratory  mechanics;  computational  model 


NEW  &  NOTEWORTHY 

Regional  ventilation  in  the  injured  lung  is  heterogeneous  and 
frequency  dependent,  making  it  difficult  to  predict  how  an 
oscillatory  flow  waveform  at  a  specified  frequency  will  be 
distributed  throughout  the  periphery.  In  this  study,  we  devel¬ 
oped  a  computational  model  of  distributed  gas  flow  and  CO2 
elimination  during  oscillatory  ventilation.  Our  simulations  in¬ 
dicate  that  ventilation  distribution  in  a  heterogeneous  lung  is 
spatially  clustered  and  frequency-dependent,  indicating  that 
oscillatory  frequency  is  an  important  factor  in  regional  gas 
exchange. 

FOR  PATIENTS  WITH  the  acute  respiratory  distress  syndrome 
(ARDS),  flow  is  distributed  throughout  the  lung  periphery  in  a 
nonuniform  and  frequency-dependent  manner  due  to  mechan¬ 
ical  heterogeneity  (3,  10).  Patients  with  ARDS  are  particularly 
susceptible  to  the  harmful  processes  of  cyclic  overdistension 
and  derecruitment,  collectively  termed  ventilator-induced  lung 
injury  (VILI).  The  ensuing  mechanical  stresses  associated  with 
these  processes  result  in  the  release  of  inflammatory  mediators 
that  may  contribute  significantly  to  multiple  organ  failure  and 
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death  (37).  Intratidal  overdistension  and  derecruitment  often 
occur  simultaneously  in  different  regions  of  injured  lung  due  to 
nonuniform  distribution  of  flow  (8).  Thus  optimizing  the  mag¬ 
nitude  and  distribution  of  flow  to  mitigate  VILI  is  a  critical 
objective  of  mechanical  ventilation  in  ARDS. 

Lung-protective  ventilation  significantly  reduces  mortality 
in  patients  with  ARDS,  using  high  positive  end-expiratory 
pressures  to  prevent  derecruitment  and  small  tidal  volumes  to 
avoid  overdistension  (1).  High-frequency  oscillatory  ventila¬ 
tion  (HFOV)  is  an  alternative  modality  that  applies  ventilatory 
rates  >2  Hz  and  is  used  as  a  rescue  therapy  for  ARDS  patients 
whose  condition  is  refractory  to  conventional  or  lung-protec¬ 
tive  ventilation  (31).  Despite  the  use  of  smaller  tidal  volumes 
compared  with  lung-protective  ventilation,  multicenter  clinical 
trials  of  HFOV  in  adult  ARDS  patients  failed  to  demonstrate 
improvements  in  mortality  (12,  56).  This  may  be  due  to 
heterogeneous  and  frequency-dependent  distributions  of  flow 
in  the  periphery,  which  may  contribute  to  regional  lung  injury 
that  is  not  easily  observed  in  patients. 

There  are  several  mechanisms  by  which  gas  transport 
occurs  during  oscillatory  flow:  direct  alveolar  ventilation, 
asymmetric  velocity  profiles,  inspiratory-expiratory  differ¬ 
ences  in  velocity  profiles,  turbulent  and  oscillatory  disper¬ 
sion,  pendelluft  phenomenon,  cardiogenic  mixing,  collateral 
ventilation,  and  molecular  diffusion  (9,  31).  The  relative 
magnitude  of  each  mechanism’s  effect  within  any  particular 
region  of  the  lung  may  depend  on  the  combination  of 
velocity  profile,  oscillatory  frequency,  and  airway  geometry 
(9,  30).  Computational  modeling  of  gas  transport  in  the  lung 
allows  for  useful  predictions  of  the  impact  of  frequency  and 
tidal  volume  on  regional  gas  transport  in  critically  ill  pa¬ 
tients.  Such  predictions  may  provide  insight  into  derange¬ 
ments  of  gas  exchange  in  mechanically  heterogeneous 
lungs.  Computational  models  exist  for  predicting  effects  of 
specific  gas  transport  mechanisms  on  total  gas  exchange 
during  either  conventional  mechanical  ventilation  (32)  or 
HFOV  (14).  Few  models  are  applicable  to  gas  transport  in  a 
branching  airway  structure  over  a  wide  range  of  frequencies 
(25,  33),  yet  these  are  limited  by  assumptions  of  symmet¬ 
rical  airway  networks  and  uniform  distributions  of  flow. 
Current  computational  models  are  limited  in  their  ability  to 
simulate  heterogeneous  gas  transport  (e.g.,  CO2  elimination) 
spanning  a  range  of  frequencies  and  tidal  volumes  encom¬ 
passing  conventional  mechanical  ventilation  and  HFOV. 

The  goal  of  this  study  was  to  develop  a  computational  model 
of  regional  flow  and  gas  exchange  in  a  three-dimensional 
airway  network  during  oscillatory  ventilation.  Specifically,  our 
objectives  were:  7)  to  generate  frequency-dependent  distribu¬ 
tions  of  advective  flow  throughout  an  anatomically  structured 
airway  network,  imposing  either  homogeneous  or  heteroge¬ 
neous  distributions  of  tissue  elastance  to  reflect  healthy  and 
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injured  parenchymal  mechanics;  2)  to  simulate  gas  exchange 
for  any  given  tidal  volume,  oscillatory  frequency,  and  flow 
distribution,  based  on  models  of  diffusive  and  advective  gas 
transport  phenomena;  and  Anally  3)  to  assess  the  frequency- 
dependence  of  tidal  volumes  required  for  eucapnic  ventilation, 
as  well  as  ventilation  heterogeneity  for  healthy  and  injured 
lungs. 
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Fig.  1 .  Central  aii^way  s  of  a  canine  lung  as  segmented  from  an  X-ray  computed 
tomographic  scan  in  a  supine  dog  {left)  and  resulting  anatomically  based 
3-dimensional  airway  network  {right)  generated  using  a  space-filling  algorithm 
(38,  39). 
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Model  structure.  The  model  structure  is  based  on  a  three- 
dimensional  network  of  a  canine  airway  tree  shown  in  Fig.  1,  in 
which  each  airway  segment  was  assumed  to  be  an  ideal  cylinder 
with  dimensions  and  spatial  position  based  on  a  thoracic  CT  scan 
of  a  supine  dog  inflated  to  total  lung  capacity  (22).  Airways  <2 
mm  diameter  were  algorithmically  generated  within  the  con¬ 
strained  boundaries  of  segmented  lung  lobes  (38).  The  length  (/) 
and  diameter  {d)  of  each  segment  were  then  scaled  down  from  total 
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lung  capacity  to  a  mean  airway  opening  pressure  of  10  cmH20, 
where  the  size  scaling  for  each  airway  was  determined  according  to 
its  time-averaged  local  transmural  pressure,  defined  as  the  differ¬ 
ence  between  mean  lumenal  pressure  and  local  pleural  pressure  (3, 
41).  Local  pleural  pressure  in  the  model  was  assumed  to  vary  with 
gravity  by  0.25  cmH20/cm  according  to  the  hydrostatic  weight  of 
pleural  fluid  in  a  supine  dog  (26),  with  an  average  value  of  -5 
cmH20.  Airway  dimensions  defined  longitudinal  and  shunt  imped¬ 
ances  to  oscillatory  flow,  which  are  briefly  summarized  here  from 
Colletti  et  al.  (10).  Airway  segment  longitudinal  impedance  was 
determined  assuming  Womersley-type  oscillatory  flow  in  a  cylin¬ 
der  (42,  54).  Airway  segment  shunt  impedance  was  given  by  a 
parallel  combination  of  isothermal  gas  compression  compliance, 
soft  tissue  viscoelasticity,  and  cartilaginous  viscoelasticity  (10, 
24).  Each  terminal  airway  in  the  tree  was  subtended  by  a  viscoelas¬ 
tic  “constant-phase”  element  representing  acinar  wall  distensibility 
(17)  in  parallel  with  gas  compression  compliance  (10,  24).  The 
arrangements  of  all  longitudinal  and  shunt  impedances  in  the 
model  were  defined  according  to  Kaczka  et  al.  (24).  Flow  oscilla¬ 
tions  at  the  airway  opening  were  distributed  throughout  the  lung 
periphery  using  a  recursive  flow-dividing  scheme,  apportioning 
flows  at  each  topological  branching  point  among  parallel  pathways 
according  to  their  respective  total  input  impedances  (10). 

Acinar  viscoelasticity  was  determined  according  to  a  second- 
order  polynomial  function  of  local  transpulmonary  pressure  (3). 
Heterogeneous  injury  was  simulated  by  modulating  acinar  tissue 
elastance  according  to  a  gravitationally  weighted  gradient  noise 
distribution  using  a  Perlin  method  (11),  which  produces  a  textured 
distribution  that  is  smoother  and  more  natural  compared  with 
purely  random  variations  in  tissue  properties  (10).  Gravitational 
weighting  increased  the  magnitude  of  the  noise  distribution  toward 
the  dependent  regions  of  the  model,  loosely  mimicking  the  effects 
of  fluid  accumulation  in  injured  lungs.  The  resulting  distribution  of 
injured  acinar  tissue  elastance  demonstrated  a  mean  ±  standard 
deviation  of  1.93  ±  1.19  cmH20/p.l,  compared  with  the  healthy 
case  of  0.31  ±  0.02  cmH20/p.l.  Cross-sectional  views  of  acinar 
tissue  elastance  are  provided  in  Fig.  2. 

To  simplify  computation  of  mechanical  impedance  and  flow 
distribution  in  the  frequency-domain  we  did  not  consider  intratidal 
variations  in  model  structure  or  mechanical  properties  (2).  Addi¬ 
tionally,  CO2  transport  throughout  the  airway  network  was  solved 
as  a  steady-state  problem,  i.e.,  time-averaged  throughout  the  os¬ 
cillatory  ventilation  cycle  and  with  CO2  elimination  equal  to  CO2 
production. 

Dijfusive  transport.  To  approximate  CO2  diffusion  in  the  tree,  we 
assumed  that  molar  flux  of  CO2  in  the  axial  direction  of  each  airway 

segment  due  only  to  diffusive  phenomena  (M‘^^C02)  was  propor¬ 
tional  to  its  axial  partial  pressure  difference  (APsegC02): 


M^COj  = 


r^PseoCOy 


where  is  the  diffusive  transport  resistance: 


R 


diff 

T 


4TZTI 


(1) 


(2) 


The  TZ  parameter  is  the  universal  gas  constant,  T  is  body  temper¬ 
ature,  and  Deff  the  effective  diffusivity  constant  for  CO2  contained  in 
the  airway  segment.  Under  zero-flow  conditions,  Deff  is  equal  to  the 
molecular  diffusivity  Dmoi-  However  during  advective  flow,  the  ef¬ 
fective  diffusivity  Deff  includes  dispersive  components  to  account  for 
enhanced  diffusion  resulting  from  transient  velocity  profile  changes 
during  laminar  oscillatory  flow  (50)  or  the  mixing  of  eddies  during 
turbulent  oscillatory  flow  (14,  40): 


r  f^rnol  “h  ^dis  ’  RScrit 

where  Djj|”  and  Djj|*  are  the  laminar  and  turbulent  oscillatory  disper¬ 
sion  coefficients,  respectively,  and  Re  is  the  Reynolds  number  given 
by  the  airway  segment  diameter  d,  kinematic  gas  viscosity  v,  and 
root-mean-square  mean-axial  velocity  Drms: 


Re 


V 


(4) 


The  critical  Reynolds  number  (Recnt  =  30)  characterizes  the 
transition  from  laminar  to  turbulent  oscillatory  flow  (35).  While  Djf” 
is  dependent  on  oscillatory  frequency  (35,  50),  D^“f  depends  only  on 
the  root-mean-square  mean-axial  velocity  that  characterizes  eddy 
mixing  strength  (14).  To  account  for  additional  advective  gas  mixing 
at  airway  bifurcations,  we  developed  a  mixing  transport  resistance 
(R™)  for  each  segment: 


R 


mix 

T 


TZT 


(5) 


where/is  the  oscillatory  frequency  of  the  advective  flow  waveform  in 
Hz,  and  Umix  is  the  amplitude  of  volume  flux  through  the  airway 
segment.  Thus,  characterizes  net  gas  transport  from  mixing  the 
volume  of  gas  passing  completely  through  an  airway  segment  during 
each  period  of  oscillatory  flow.  Therefore  7?™“  will  be  nonzero  only 
when  this  oscillating  volume  exceeds  the  volume  of  the  airway 
segment,  Vseg,  given  by: 


V  >V 

*  mix  *  seg 


4 


where  Vmix  is: 


l^seg(/)l 


-  V 


(6) 


(7) 


Acinar  tissue  stiffness  [cmH20  pL'] 


0.26  0.58  1.26  2.74  5.98 


I  Gravity 
I  direction 


Apex  ^  ^  Base 

Fig.  2.  Transverse  cross  sections  of  simulated  regional  acinar  tissue  stiffness 
for  healthy  (upper)  and  injured  (lower)  conditions,  spatially  interpolated 
between  acinar  elements  for  visualization. 
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and  IVseg(/)l  is  the  magnitude  of  oscillatory  flow  in  the  airway  segment 
at  frequency  /.  The  mixing  transport  resistance  is  derived  in  more 
detail  in  appendix  A-1. 

The  total  transport  resistance  for  each  airway  (i?“)  is  the  parallel 
combination  of  diffusive  and  mixing  transport  resistances: 


R 


tot 

T 


1  ndiff  Dinix 

1  /V-p 

r  “ 

odiff  r>mix 

/Vp  JVp 


(S) 


An  upper  bound  on  is  obtained  in  the  zero-flow  limiting  case  when 
Vmix  ^  0  and  1?““  ^  In  this  case  =  1?“  and  Oeff  =  Omoi- 
Once  is  computed  for  each  airway  in  the  network,  a  mesh 
analysis  (appendix  A-2)  forms  a  system  of  equations  for  solving  the 
molar  flux  of  CO2  at  each  acinus,  given  the  partial  pressure  differen¬ 
tial  between  each  acinus  and  the  airway  opening: 


R^MCOi  =  APCO2  (9) 


where  R-r  is  a  full  matrix  of  transport  resistances  (appendix  A-2), 
MCO2  is  a  vector  of  acinar  CO2  molar  fluxes  (i.e.,  for  the  nth  acinus: 
M„C02),  and  APCO2  is  a  vector  of  CO2  partial  pressure  differentials 
between  acini  and  airway  opening  (i.e.,  for  the  nth  acinus:  AP„C02  = 
P„C02  —  /'aoC02).  Each  row  of  the  system  corresponds  to  the  mesh 
formed  by  the  airways  connecting  a  single  acinus  to  the  airway 
opening.  The  partial  pressure  of  CO2  at  the  airway  opening  was 
assumed  to  be  zero.  After  solving  the  system  of  equations,  the  CO2 
elimination  due  to  diffusion,  dispersion,  and  mixing  (V'™'^C02)  for  each 
acinus  can  be  converted  from  the  molar  flux  using  the  ideal  gas  law: 

TZT 

^yODMco^  =  —M„C02  {10) 

atm 


The  mesh  analysis  ensures  conservation  of  mass  between  CO2  elim¬ 
ination  at  the  airway  opening  (V™^C02)  and  the  sum  of  V°™C02 
across  all  acini. 

Advective  transport.  The  distribution  of  advective  flow  throughout 
the  airway  tree  was  computed  using  a  recursive  flow  divider  algorithm 
(10).  Direct  acinar  ventilation  occurs  when  volume  amplitudes  at  the 
airway  opening  are  sufficiently  large  to  allow  fresh  gas  mixing  directly 
with  acinar  gas  during  inspiration,  effectively  bypassing  and  in 
Eq.  8.  To  model  this  mechanism  of  gas  transport,  we  assume  that  the 
volume  of  fresh  gas  delivered  to  each  acinus  is  completely  mixed  during 
inspiration  and  that  the  same  volume  of  mixed  gas  is  completely  removed 
during  expiration,  according  to  Fig.  3.  CO2  elimination  due  to  direct 
acinar  ventilation  is  computed  according  to  the  concept  of  “personal  dead 
space,”  which  apportions  an  amount  of  inspired  fresh  gas  to  each 
individual  acinus  (13). 

In  this  manner,  the  dead  space  of  all  conducting  airways  was  succes¬ 
sively  distiibuted  in  proportion  to  the  flow  distribution  thi'ough  subtend¬ 
ing  airways  until  reaching  the  terminal  acini.  Thus  the  distribution  of 
personal  dead  space  is  dependent  on  frequency,  in  accordance  with  the 
distribution  of  flow.  We  determined  the  CO2  elimination  due  to  direct 
acinar  ventilation  (V^^''C02)  for  a  particular  acinar  compartment  by 
comparing  its  personal  dead  space  volume  [VD,n(/)]  to  the  total  volume 
delivered  to  that  acinus  over  one  period  of  oscillation: 

P  CO 

v;?^''co2  =  ^/|vp,„|  (11) 


where  863  is  a  conversion  factor  between  partial  pressure  in  mmHg  and 
gas  volume  fraction,^  PnCOi  is  the  partial  pressure  of  acinar  CO2,  and 


^  The  conversion  factor  accounts  for  the  difference  between  standard  tem¬ 
perature  and  pressure  dry  (STPD)  air  and  body  temperature  and  pressure  air 
saturated  with  water  vapor  (BTPS). 


End-expiration 


End-inspiration 


Fig.  3.  Direct  acinar  ventilation  is  modeled  by  assuming  complete  mixing 
within  each  nth  terminal  acinus  of  inspired  fresh  gas  (FF,n),  inspired  dead  space 
gas  (Vd,/j),  and  residual  volume  (Va,/i)-  The  concept  of  “personal  dead  space” 
is  used  to  distribute  the  total  dead  space  (Vd)  throughout  all  acini  in  proportion 
to  the  magnitude  of  the  ventilation  distribution  (13). 


IVf,«I  is  the  magnitude  of  frequency-dependent  fresh  gas  volume  deliv¬ 
ered  to  the  acinus  given  by: 


^F,.|  - 


Uf)\ 


-  \VD.nif)\, 


\Uf)\ 


> 


I  ^D,„(/)  I 


(12) 


0, 


otherwise 


\V„{f)\  is  the  magnitude  of  oscillatory  flow  at  frequency  /  in  the 
terminal  acinus.  Note  that  if  the  total  volume  of  inspired  gas  delivered 
to  an  acinus  over  one  period  of  oscillation  is  less  than  the  personal 
dead  space  of  that  acinus,  there  will  be  no  contribution  of  direct  acinar 
ventilation  to  CO2  elimination  from  that  particular  acinus.  Total  dead 
space  volume  (Vd)  was  calculated  as  the  sum  of  all  cylindrical  airway 
segment  volumes. 

Total  transport  and  eucapnia.  Total  CO2  elimination  (V'totC02)  for 
the  model  at  a  specified  frequency  was  computed  from  the  sum  of 
CO2  elimination  across  all  N  acini: 


VtotC02  =  2  %C02 


(13) 


where  n  indexes  individual  acini,  and  the  CO2  elimination  in  each 
acinus  (V'„C02)  is  determined  by  the  sum  of  CO2  elimination  via 
direct  acinar  ventilation  and  by  solving  the  system  of  transport 
resistances: 


1/„C02  =  V^^''C02  +  \>°™C02 


PnCO: 

863 


■f\yF.n\ 


TZT  . 

+  — M„C02 

^atm 


(14) 


Equation  14  assumes  that  each  acinar  CO2  partial  pressure  is 
equal  to  its  corresponding  end-capillary  CO2  partial  pressure. 
Mixed-venous  blood  was  assumed  to  enter  each  capillary  with  46 
mmlTg  CO2  partial  pressure.  However  end-capillary  CO2  partial 
pressure  for  each  acinus  is  determined  by  the  rates  of  acinar  CO2 
elimination  and  perfusion.  Thus  the  distribution  of  P„C02  is 
dependent  on  the  distribution  of  acinar  CO2  elimination,  and  Eq. 
14  for  total  CO2  elimination  must  be  solved  iteratively.  An  opti¬ 
mization  routine  minimized  the  error  in  conservation  of  mass  for 
each  acinar  compartment,  such  that  CO2  elimination  via  ventilation 
(i.e.,  Eq.  14)  was  equal  to  the  product  of  perfusion  rate  and  CO2 
content  difference  between  end-capillary  and  mixed-venous  blood 
(53).  Eucapnic  ventilation  was  assumed  to  occur  when  total  CO2 
elimination  for  the  model  approximately  equaled  the  predicted 
metabolic  CO2  production  for  a  25  kg  dog  (1.93- 10^^-Ls^')  (18). 

Simulations.  The  recursive  algorithms  for  computing  and  storing 
impedances  and  determining  advective  flow  division  throughout  the 
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Tidal  volume  [  L  ] 

Fig.  4.  Simulation  results  for  ytotC02  as  a  function  of  Vt  and/.  Solid  black 
horizontal  line  represent  the  lower  bound  of  total  CO2  elimination  due 
solely  to  molecular  diffusion  (i.e.,  zero-flow  conditions),  while  the  dashed 
horizontal  line  represents  total  CO2  elimination  to  achieve  eucapnia  for  a 
25  kg  dog.  Text  labels  indicate  dominant  gas  transport  mechanism  during 
0.1  Hz  oscillation. 

tree  were  written  in  C  +  + .  Computation  time  for  generating  advective 
flow  distribution  at  each  distinct/was  ~3  s  on  a  64-bit  computer  with 
an  Intel  Core  i7-950  processor  operating  at  3.07  GHz  with  12  GB 
RAM.  The  gas  transport  model  was  written  and  executed  in  MAT- 
LAB  (version  7.13,  The  Mathworks,  Natick,  MA).  Total  computation 
time  for  each  distinct /and  Vt  simulation  was  -300  s. 

The  mechanical  impedance  and  ventilation  distribution  for  the 
model  were  computed  at  115  distinct  frequencies,  spanning  0.1-100 
Hz.  The  resonant  frequency  (/-es),  at  which  the  reactive  component  of 
total  lung  impedance  is  zero,  was  determined  hy  interpolation.  Gas 
transport  distributions  were  computed  for  a  range  of  volume  ampli¬ 
tudes  at  each  frequency  between  0.1  and  30  Hz,  with  emphasis  on 
smaller  volume  amplitudes  at  higher  frequencies.  Volume  amplitudes 


were  selected  between  5  and  600  ml,  chosen  to  ensure  that  a  simula¬ 
tion  representative  of  eucapnic  ventilation  was  computed  at  each 
frequency.  The  tidal  volume  required  for  eucapnic  ventilation  (V^“)  as 
a  function  of  frequency  was  regressed  to  a  power-law  (47): 

VTif)  =  Po  +  i.15) 

The  parameters  (So,  Pi,  and  P2  were  estimated  using  a  nonlinear 
least-squares  technique  (MATLAB  v7.13). 

To  quantify  the  heterogeneity  of  ventilation  distribution,  we  exam¬ 
ined  the  amplitude  and  phase  histograms  of  acinar  flows  as  functions 
of  oscillatory  frequency  for  both  health  and  injured  conditions,  along 
with  the  histograms  for  acinar  CO2  elimination.  Acinar  flow  magni¬ 
tudes  were  normalized  relative  to  the  theoretical  values  obtained  for  a 
perfectly  symmetric,  homogeneous  lung  with  rigid  airway  walls  and 
no  gas  compression  (i.e.,  perfectly  uniform  ventilation  distrihution). 
Acinar  flow  phases  were  normalized  relative  to  the  tracheal  flow 
phase,  constrained  to  be  within  ±  180°.  Bin  sizing  for  normalized 
acinar  flow  magnitudes,  phases,  and  V'„C02  were  167  per  decade,  1.39 
per  degree,  and  167  per  decade,  respectively. 

RESULTS 

Validation.  Figure  4  shows  simulation  results  for  ytotC02  as 
a  function  of  Vt  and  /  during  oscillatory  ventilation  with  a 
single  frequency.  Simulations  are  shown  for  each  frequency 
with  increasing  V/  increments  of  5  ml,  until  14^,^002  exceeded 
the  predicted  value  for  eucapnia.  The  required  CO2  elimination 

to  maintain  eucapnia  (y(“(‘^C02)  is  also  indicated,  as  well  as 
Vt(,tC02  produced  by  molecular  diffusion  only  (i.e.,  zero-flow 
conditions)  which  establishes  the  lower  bound  of  CO2  elimi¬ 
nation  regardless  of  Vt  or  /.  There  are  noticeably  different 
behaviors  for  VjjjjC02  with  respect  to  Vt  and  /.  Specifically, 
there  are  certain  ranges  of  volume  amplitudes  for  which  certain 
gas  transport  mechanisms  dominate  the  others:  starting  from 
the  lower  bound  molecular  diffusion,  there  are  dramatic  in¬ 
creases  in  Vjq,C02  with  the  advent  of  oscillatory  dispersion, 
followed  by  turbulent  dispersion,  advective  mixing  at  bifurca¬ 
tions,  and  finally  direct  acinar  ventilation.  These  ranges  were 
identified  by  iteratively  eliminating  specific  gas  transport 
mechanisms  from  the  model,  and  then  resimulating  VtotC02. 

Figure  5  shows  as  a  function  of  frequency  over  local¬ 
ized  frequency  domains  0. 1  Hz  <  /  <  1  Hz  and  2  Hz  <  /  <  30 
Hz.  Regression  curves  are  extended  beyond  the  fit  domains  to 
demonstrate  the  relative  reductions  in  required  volume  ampli¬ 
tudes  that  occur  inside  each  range  compared  with  the  extrap¬ 
olated  fit  of  the  other  range.  Parameter  estimates  are  provided 
in  Tables  1  and  2.  Expected  values  for  low-frequency  gas 
exchange  (i.e.,  0.1  Hz  </<  1  Hz)  are  derived  from  standard 
equations  of  gas  exchange  during  conventional  mechanical 


Fig.  5.  Simulation  results  representing  model 
predictions  for  eucapnic  ventilation  at  vari¬ 
ous  single  frequencies  of  oscillation  (dashed 
black  lines).  Dotted  vertical  black  lines  rep¬ 
resent  the  resonant  frequency  (fx&s-  elastic 
and  inertial  components  of  impedance  equal 
and  opposite).  Model  predictions  of  tidal 
volume  required  for  eucapnic  ventilation  de¬ 
rived  from  power-law  regressions  over  low 
frequencies  (purple:  0.1  ^  ^  1  Hz)  and 

high  frequencies  (blue:  2  </<  30  Hz). 
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Table  1.  Power-law  regression  parameters  for  Eq.  15 
during  low  frequency  oscillatory  ventilation 


Parameter 

Expected  Value 

Regression 

Estimate 

Deviation  from 
Expected  Value,  % 

Po 

Vo  =  0.169 

0.151 

10.7 

Pi 

863  . 

- Vf/C02  =  0.042 

PF-Oi 

0.043 

2.2 

P2 

-1 

-0.99 

1.1 

ventilation  and  spontaneous  breathing  (appendix  A-3).  Ex¬ 
pected  values  for  high  frequency  gas  exchange  (i.e.,  2  Hz  < 
/  <  30  Hz)  are  obtained  from  previous  experimental  studies. 

Figure  6  shows  a  comparison  between  these  interpolated 
eucapnic  simulation  results  and  experimental  data  from  Ven¬ 
egas  et  al.  (47),  converted  to  dimensionless  equations  of  gas 
exchange: 


/vr(/) 

{16) 

Va 

fVu 

{17) 

Va 

where  Q  and  T  are  the  normalized  flow  cost  and  normalized 
frequency,  respectively,  and  Va  is  the  total  acinar  ventilation 
required  for  eucapnia: 


'^A  = 


vzc^\ 

PiC02  I 


(18) 


where  P3CO2  is  the  eucapnic  partial  pressure  of  CO2  in  mixed 
arterial  blood.  The  low-frequency  relationship  Q  =  \  -\-  J- 
corresponds  to  theoretical  gas  exchange  primarily  via  direct 
acinar  ventilation,  whereas  the  high-frequency  relationship 
Q  =  is  a  regression  to  experimentally  acquired  data 

(47).  The  model  accurately  predicts  the  transition  from  low- 
frequency  to  high-frequency  gas  exchange  behavior,  and  sim¬ 
ulation  results  demonstrate  close  qualitative  and  quantitative 
agreement  with  the  experimental  data  of  Venegas  et  al.  (47), 
deviating  at  the  highest  frequencies  {IF  >  60,  or/>  12.8  Hz). 

Impact  of  heterogeneity.  Figure  7  shows  five  transverse 
model  cross  sections  that  are  color-coded  to  depict  the  distri- 


Table  2.  Ratio  of  power-law  exponents  for  CO2  elimination 
as  a  function  off"  and  V*  during  high  frequency  oscillatory 
ventilation,  as  reported  in  various  studies 


Study 

Species 

h/a 

Slutsky  et  al.  1981  (36)* 

dog 

1.7 

Berdine  et  al.  1984  (7)* 

dog 

1.9 

Jaeger  1991  (19) 

dog 

1.9 

Jaeger  et  al.  1984  (20) 

dog 

2.0 

Venegas  et  al.  1985  (45) 

dog 

2.1 

Present  study 

dog 

2.1 

Venegas  et  al.  1986  (47) 

various 

2.2 

Watson,  Jackson  1985  (52) 

monkey 

1.9 

Watson,  Jackson  1984  (51) 

rabbit 

2.7 

Only  the  ratio  of  b  to  a  is  presented  because  the  power-law  regression  in  Eq. 
J5  determines  only  the  relative  relationship  between  frequency  and  volume 
amplitude.  *Values  for  a  and  b  reported  by  Venegas  et  al.  (47).  fValues  for  a 
and  b  reported  by  Chang  (9). 


Fig.  6.  Comparison  between  simulation  results  and  experimental  data  of 
Venegas  et  al.  (47).  The  axes  are  dimensionless  variables  of  frequency  (F)  and 
flow  (Q).  The  black  lines  correspond  to  theoretical  predictions  of  low- 
frequency  gas  exchange  (dashed)  and  regression  to  experimental  HFOV  data 
(solid).  The  gray  line  represents  the  simulations  results  interpolated  at  eucapnic 
conditions. 


bution  of  normalized  acinar  flows,  along  with  acinar  CO2 
elimination  rates,  for  frequencies  between  0.2  and  30  Hz. 
Acinar  flow  magnitudes  were  normalized  relative  to  their 
corresponding  values  for  a  perfectly  symmetric,  homoge¬ 
neous  lung  with  rigid  airway  walls  and  no  gas  compression. 
The  distribution  of  acinar  CO2  elimination  was  normalized 
to  the  total  CO2  elimination  divided  by  the  number  of  acini. 
The  coefficients  of  determination  between  acinar  flow  mag¬ 
nitude  and  CO2  elimination  (r^)  at  each  frequency  and  condi¬ 
tion  are  also  shown.  Acinar  CO2  elimination  was  positively 
correlated  with  acinar  flow  magnitude  for  either  condition, 
although  the  correlation  tended  to  decrease  with  increasing 
frequency.  Gravitational  gradients  were  apparent  at  low  fre¬ 
quencies  for  the  healthy  condition,  yet  diminish  with  increas¬ 
ing  frequency.  For  the  injured  condition,  the  distributions  of 
both  flow  and  CO2  elimination  were  spatially  concordant  with 
the  distribution  of  acinar  tissue  stiffness  (Fig.  2)  for  the 
frequencies  shown. 

Figure  8  shows  two-dimensional  histograms  of  normalized 
acinar  flow  magnitude  and  phase  distributions  for  healthy  and 
injured  conditions,  along  with  the  acinar  CO2  elimination 
distribution  from  0.1  to  100  Hz.  Acinar  flow  phase  was 
normalized  relative  to  the  corresponding  tracheal  flow  phase, 
constrained  to  be  within  ±  180°.  The  histograms  are  color- 
coded  based  on  the  percentage  of  acini  in  the  model.  For 
frequencies  less  than  resonance,  ventilation  distribution  is 
largely  frequency  independent.  However,  when  the  oscillatory 
frequency  exceeds  resonance,  ventilation  distribution  becomes 
both  frequency  dependent  and  heterogeneous,  even  for  the 
mechanically  homogeneous  healthy  case.  The  low-frequency 
ventilation  distribution  is  much  more  homogeneous  for  the 
healthy  case  compared  with  the  injured  case,  as  expected. 
However  the  effect  of  high  frequency  on  ventilation  heteroge¬ 
neity  is  comparable  between  the  two  conditions,  in  both  cases 
increasing  heterogeneity  with  increasing/ > /res.  Above /es  the 
mean  normalized  acinar  flow  magnitude  exceeded  unity,  indi- 
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Fig.  7.  Transverse  cross  sections  of  simulated  ^ 
regional  acinar  flow  magnitude  {left)  and  CO2 
elimination  {right)  at  selected  oscillatory  fre¬ 
quencies,  for  healthy  {top)  and  injured  {bottom) 
conditions,  spatially  interpolated  between  aci¬ 
nar  elements  for  visualization.  Coefficients  of 
determination  (r^)  between  acinai*  flow  magni¬ 
tude  and  CO2  elimination  are  shown  at  the  far 
right. 
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eating  that  the  sum  of  acinar  flow  magnitudes  was  greater  than 
the  tracheal  flow  magnitude.  The  increase  in  mean  acinar  flow 
magnitudes  above /res  is  coincident  with  increased  acinar  phase 
variance.  Multimodal  patterns  were  also  visible  in  histograms, 
reflecting  large  groups  of  acini  with  similar  frequency-depen- 
dent  behavior. 

To  determine  whether  such  multimodality  in  ventilation 
distribution  was  consistent  with  anatomic  grouping,  we  color- 
coded  individual  acini  according  to  the  frequency  at  which  they 
experienced  either  a  local  minimum  (fmin)  or  maximum  (/max) 
normalized  flow  magnitude.  Figure  9  shows  the  spatial  ana¬ 
tomic  distributions  of  /min  and  /max  for  the  three-dimensional 
model.  Spatially  organized  clusters  of  acini  with  the  same /min 
or /max  values  are  immediately  apparent,  implying  that  over¬ 
ventilation  and  underventilation  may  manifest  simultaneously 
in  clusters  of  the  lung,  in  a  frequency-dependent  manner.  The 
particular  pattern  of  clustering  was  similar  between  the  healthy 
and  injured  conditions,  despite  substantial  differences  in  the 
distribution  of  tissue  elastance  (Fig.  2). 

DISCUSSION 

In  this  study,  we  developed  an  anatomic  computational 
model  of  gas  exchange  during  oscillatory  ventilation  by  incor¬ 
porating  mechanical  impedance  and  pertinent  gas  transport 
mechanisms  over  a  wide  range  of  frequencies.  The  dominant 
mechanism  of  gas  transport  exhibits  strong  scaling  dependence 
on  Vt,  with  transitions  from  molecular  diffusion  near  zero-flow 
conditions  to  laminar  oscillatory  dispersion,  then  turbulent 
dispersion,  then  advective  mixing  at  bifurcations,  and  finally 
direct  acinar  ventilation  at  the  largest  V/-  The  transition  from 


laminar  dispersion  to  turbulent  dispersion  is  determined  by  a 
critical  Reynolds  number,  which  was  set  at  Recut  =  30  (35). 
Values  of  Re  are  proportional  to  the  root-mean-square  axial 
flow  velocity,  which  implies  dependence  on  both  volume 
amplitude  and  oscillatory  frequency.  Therefore,  the  transition 
from  laminar  to  turbulent  dispersion  will  occur  at  different 
volume  amplitudes  for  each  frequency.  Laminar  oscillatory 
dispersion  alone  is  not  sufficient  to  provide  eucapnic  ventila¬ 
tion  at  a  reasonably  small  V/  for  any  /  in  the  simulated  range 
between  0.1  and  30  Hz.  By  contrast,  the  transition  to  turbu¬ 
lence  yields  about  an  order-of-magnitude  increase  in  ytotC02, 
such  that  for  high  frequencies  (f>  10  Hz),  adequate  ytotC02  is 
achieved  primarily  via  dispersion. 

For  intermediate  frequencies  (between  1  and  10  Hz),  mixing 
at  bifurcations  contributes  the  most  to  eucapnic  ventilation, 
nearly  doubling  ytotC02  compared  with  dispersion  alone.  The 
transition  from  dispersion-dominated  to  bifurcation-mixing- 
dominated  gas  transport  occurs  at  a  constant  V/  between  50 
and  60  ml  across  all  frequencies.  The  trachea  of  our  canine 
model  has  a  lumenal  volume  of  56.3  ml,  which  is  by  far  the 
largest  volume  of  any  airway  in  the  tree:  the  next  largest  is  6.5 
ml,  followed  by  2.2  ml.  Thus  we  may  safely  assume  that 
advective  mixing  at  airway  bifurcations  in  the  model  has  the 
greatest  effect  on  total  gas  transport  for  V/  larger  than  tracheal 
volume. 

The  final  transition  in  dominant  gas  transport  mechanisms  is 
attributable  to  the  increasing  prevalence  of  direct  acinar  ven¬ 
tilation.  Occurring  only  at  the  largest  V/,  this  mechanism 
requires  a  front  of  fresh  gas  to  move  directly  between  the 
airway  opening  and  the  acinar  compartments,  producing  sub- 
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Fig.  8.  Simulation  representing  model  predictions  for  eucapnic  ventilation  vs.  oscillatory  frequency.  Horizontal  axes  are  aligned,  dashed  black  lines  represent  the 
resonant  frequency  (/ies:  elastic  and  inertial  components  of  impedance  equal  and  opposite).  Acinar  flow  magnitude,  flow  phase,  and  CO2  elimination  are 
normalized  as  described  in  the  text. 


stantial  gas  transport  via  advection.  The  transition  to  direct 
acinar  ventilation  occurs  over  a  wider  range  of  VV  compared 
with  the  other  mechanisms,  due  to  differences  in  the  ventilation 
distribution  as  well  as  the  distribution  of  path  lengths.  The 
impact  of  direct  acinar  ventilation  on  gas  transport  is  also  much 
larger  than  that  of  any  other  mechanism  for  CO2  elimination  by 
almost  an  order  of  magnitude. 

The  last  two  transitions  seem  to  have  constant  Vr  thresholds, 
whereas  the  diffusion  and  dispersion  transitions  are /dependent  as 
well  as  Vt  dependent.  This  can  be  explained  by  the  physical 
dependence  of  direct  acinar  ventilation  and  bifurcation  mixing  on 
the  size  and  arrangement  of  dead  space  in  the  conducting  airways, 
whereas  dispersion  depends  more  strongly  on  the  magnitude, 
frequency,  and  cross-sectional  profile  of  the  flow  waveform. 

Pendelluft  occurs  when  lung  regions  with  differing  mechan¬ 
ical  time  constants  oscillate  out  of  phase  with  respect  to  each 
other,  resulting  in  increased  intrapulmonary  gas  mixing  and 
homogenization  (16,  29).  In  this  model  of  time-averaged  gas 
transport,  pendelluft  manifests  as  increased  acinar  phase 
variance  and  flow  magnitude,  as  shown  in  Fig.  8.  Gas 
mixing  is  increased  between  lung  regions  oscillating  out  of 
phase,  due  to  increased  flow  magnitude  in  each  region.  The 
simulations  performed  in  this  study  indicate  that  pendelluft 
occurs  primarily  above  the  resonant  frequency,  consistent 
with  experimental  findings  of  phase  variance  using  strobo¬ 


scopic  measurements  of  the  lung  surface  (27)  and  alveolar 
pressure  capsules  (15). 

Validation.  Simulations  from  this  canine  lung  model,  under 
healthy  conditions  with  near  homogeneous  mechanical  prop¬ 
erties,  exhibit  agreement  with  experimental  results  regarding 
the  behavior  of  ytotC02  with  respect  to  /  Vr,  and  Vb  (47). 
Low-/  gas  exchange  is  dominated  by  direct  acinar  ventilation. 
The  personal  dead  space  and  alveolar  mixing  models  demon¬ 
strate  concordance  with  expected  behavior  of  CO2  elimination 
during  conventional  mechanical  ventilation  and  spontaneous 
breathing.  Theoretical  predictions  and  empirical  findings  for 
power-law  behavior  of  high-/  CO2  elimination  are  relatively 
more  inconsistent  by  comparison.  In  general  it  is  accepted  that 
VtotC02  depends  more  strongly  on  V/  than  on /according  to: 

y,otC02=^ry^  (79) 

with  a  <  b  (31).  The  exponent  values  obtained  from  these 
simulations  are  within  the  range  of  values  provided  in  the 
literature,  represented  in  Table  2  (9,  47). 

Considering  the  normalized  representation  of  flow  cost  vs. 
frequency  used  by  Venegas  et  al.  (47),  the  simulation  results 
are  in  close  qualitative  and  quantitative  agreement  with  exper¬ 
imental  results.  However,  our  simulations  underestimated  the 
normalized  flow  cost  Q  at  the  highest  frequencies  (i.e.,  T  >  60, 
or/  >  12.8  Hz)  in  animals.  This  discrepancy  may  be  attributed 
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Fig.  9.  A:  frequencies  of  maximum  and  min¬ 
imum  normalized  flow  magnitude,  shown 
for  clarity  in  the  context  of  a  single  acinus. 
B:  frequencies  of  maximum  (top)  and  mini¬ 
mum  (bottom)  normalized  flow  magnitude, 
represented  spatially  for  every  acinus,  using 
consistent  coloring  according  to  the  scale 
provided  in  A.  Acinar  flow  magnitude  is 
normalized  to  the  uniform  distribution  (unity 
corresponds  to  tracheal  magnitude  divided 
by  total  number  of  acini).  Acini  that  do  not 
exhibit  a  local  minimum  or  maximum  within 
the  range  of  frequencies  simulated  are 
shown  in  black.  The  gravitational  vector 
points  out  of  the  page,  toward  the  reader. 
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to  an  overemphasized  effect  of  turbulent  dispersion  in  the 
model  and  may  be  improved  by  the  fine-tuning  of  a  scaling 
parameter  used  to  calculate  DJjI*.  In  previous  studies,  Fredberg 
(14)  used  an  order-of-magnitude  approximation  to  set  the 
scaling  parameter  equal  to  unity,  while  Slutsky  et  al.  (35)  used 
an  average  value  obtained  from  inspiratory  and  expiratory 
flows  in  a  mechanical  model  (34)  used  an  average  value 
obtained  from  inspiratory  and  expiratory  flows  in  a  mechanical 
model.  It  is  possible  that  more  accurate  values  of  the  scaling 
parameter  during  the  turbulent  flows  of  HFOV  may  be  ob¬ 
tained  either  through  experiment  or  through  computational 
fluid  dynamics  (CFD)  simulation. 

Impact  of  heterogeneity  and  frequency  selection  during 
HFOV.  We  can  make  several  predictions  about  the  impact  of 
frequency  selection  on  ventilation  heterogeneity  during  HFOV, 
by  jointly  observing  the  frequency  dependence  of  the  acinar 
flow  distribution,  acinar  y„C02  distribution,  and  total  eucapnic 
volume  amplitude.  The  resonant  frequency  /res  marks  the 
transition  from  elastically  dominated  to  inertially  dominated 
flow.  For  frequencies  less  than /res,  flow  distribution  heteroge¬ 
neity  is  frequency  independent  and  determined  primarily  by  the 
distribution  of  resistive  and  elastic  properties.  For  frequencies 
greater  than  /es,  the  transition  to  frequency-dependent  flow 
heterogeneity  suggests  that  local  transitions  from  elastically 
dominated  to  inertially  dominated  flow  occur  at  different  fre¬ 
quencies  throughout  the  lung,  depending  on  local  mechanical, 
geometrical,  and  topological  properties.  The  multimodal  be¬ 
havior  of  the  flow  heterogeneity  (Fig.  8)  and  the  spatial 
distributions  of  /„!„  and /max  (Fig.  9)  suggest  that  frequency- 
dependent  behavior  of  the  flow  distribution  is  regionally  clus¬ 
tered.  Note  that  two  acini  sharing  the  same  value  of /max  does 
not  imply  also  sharing  the  same  value  of /min,  suggesting  that 
subtle  differences  between  mechanical  properties  of  neighbor¬ 
ing  acini  produce  noticeable  differences  in  frequency  depen¬ 
dence  of  ventilation  distribution.  Furthermore,  the  healthy  and 


injured  conditions  exhibited  similar  clustering,  suggesting  that 
the  frequency-dependent  behavior  of  ventilation  distribution  is 
largely  determined  by  airway  properties  affecting  inertance 
(i.e.,  airway  size,  path-length  between  airway  opening  and 
acinus,  airway  network  topology)  rather  than  tissue  properties 
in  this  particular  model  of  injury. 

The  transition  from  low-/  to  high-/  gas  exchange,  which 
occurs  between  0.5  and  5  Hz  in  our  canine  model,  is  associated 
with  increased  acinar  I4CO2  heterogeneity  despite  no  change 
in  flow  distribution  heterogeneity.  This  can  be  explained  by  the 
nonuniform  distribution  of  personal  dead  space  across  the 
acini.  Even  if  flow  is  distributed  uniformly,  differences  in 
path-length,  airway  network  symmetry,  and  airway  sizes  cause 
some  acini  to  receive  less  fresh  gas  per  oscillatory  period  than 
others.  As  the  rate  of  ventilation  increases,  the  tidal  volume 
required  for  eucapnia  decreases,  and  the  acini  with  large 
personal  dead  space  do  not  participate  in  direct  acinar  ventila¬ 
tion.  A  bimodal  distribution  of  acinar  CO2  elimination  exists 
over  the  transition  range  0.5  Hz  <  /  <  5  Hz,  with  one  mode 
corresponding  to  acini  receiving  direct  acinar  ventilation  and 
the  other  mode  correspond  to  acini  eliminating  CO2  only  via 
mixing  and  dispersion.  For  example,  the  fraction  of  acini 
receiving  direct  acinar  ventilation  in  the  healthy  simulations 
reduced  to  50%  at  1.6  Hz  and  0%  at  5  Hz.  In  this  range  of 
frequencies,  increased  acinar  y„C02  heterogeneity  gives  the 
appearance  of  increased  ventilation-to-perfusion  mismatching. 
The  effects  of  frequency-dependent  transitioning  from  direct 
acinar  ventilation  to  mixing  and  dispersion  observed  in  these 
simulations  are  consistent  with  the  results  of  washout  imaging 
reported  by  Venegas  et  al.  (49),  whose  findings  were  supported 
by  a  simplified  model  with  heterogeneous  dead  space.  Achiev¬ 
ing  increased  direct  acinar  ventilation  at  higher  frequencies, 
and  thereby  reducing  the  tidal  volume  required  for  eucapnia, 
may  be  possible  through  the  use  of  specialized  ventilation 
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apparatus  designed  to  minimize  dead  space  of  the  ventilator 
tubing  and  upper  respiratory  tract  (44). 

The  selection  of  an  optimal  frequency  for  delivering  HFOV 
may  be  described  as  a  multiobjective  optimization  problem, 
wherein  one  objective  is  to  minimize  the  volume  amplitude  and 
another  objective  is  to  minimize  ventilation  heterogeneity. 
Lower  frequencies  (f<fies)  require  greater  volume  amplitudes, 
thus  increasing  the  risk  of  overdistension  in  the  presence  of 
lung  injury.  High  frequencies  (f  >  /res)  minimize  the  volume 
amplitude  required  for  eucapnia  yet  increase  acinar  flow  het¬ 
erogeneity  and  consequently  acinar  1^,1002  heterogeneity. 
Based  on  these  simulations,  we  predict  that  optimal  HFOV  for 
this  canine  model  is  delivered  using  frequencies  slightly  less 
than  /res-  This  range  provides  the  best  compromise  between 
minimizing  both  volume  amplitude  and  ventilation  heteroge¬ 
neity  in  the  simulated  canine  lungs. 

The  spatial  distributions  shown  in  Fig.  7  demonstrate  spatial 
concordance  between  regional  acinar  flow  and  CO2  elimina¬ 
tion,  consistent  with  experimental  studies  combining  positron 
emission  and  X-ray  computed  tomographies  (21,  43,  55).  The 
positive  correlation  between  flow  and  CO2  elimination  de¬ 
creased  at  higher  frequencies,  likely  due  to  the  diminishing 
influence  of  direct  acinar  ventilation  and  bulk  transport  on  CO2 
elimination.  Simulations  from  the  healthy  condition  demon¬ 
strate  a  transition  from  elastically  dominated  flow  and  CO2 
elimination  at  low  frequencies  (as  indicated  by  the  influence  of 
gravity  on  these  distributions)  to  variable  regional  ventilation 
at  higher  frequencies.  That  is,  the  same  lung  region  can 
experience  a  larger  share  of  total  ventilation  at  one  particular 
frequency,  yet  smaller  shares  at  others.  This  phenomenon 
yields  the  clustering  of  regionally  preferred  frequencies /max  as 
shown  in  Fig.  9.  This  finding  suggests  that  a  combination  of 
multiple  frequencies  delivered  simultaneously  may  produce 
more  uniform  overall  ventilation  (23). 

Limitations 

Despite  the  complexity  of  our  model  in  terms  of  structural 
geometry,  topology,  fluid  flow,  and  gas  transport,  there  are 
numerous  assumptions  and  simplifications  that  limit  the  phys¬ 
iological  interpretation  of  these  results.  The  pathophysiological 
features  of  ARDS  are  not  realistically  described  by  a  simple 
altered  distribution  of  tissue  elastances.  The  model  also  does 
not  account  for  derecruitment,  perfusion  shunting,  or  any 
changes  in  airway  mechanics.  Furthermore,  the  model  ignores 
intratidal  variations  in  mechanics,  which  are  especially  rele¬ 
vant  to  overdistension  and  cyclic  derecruitment  (2,  5). 

Another  assumption,  that  mixed-venous  blood  has  a  constant 
partial  pressure  of  CO2  in  all  simulations,  restricts  the  appli¬ 
cability  of  these  simulations  to  eucapnic  conditions  and  short 
time  scales;  that  is,  all  simulations  are  executed  to  calculate  an 
instantaneous  VtotC02  under  the  premise  that  mixed-venous 
CO2  partial  pressure  is  uniformly  46  mmHg  throughout  the 
lung.  If  the  value  of  instantaneous  VtotC02  calculated  in  this 
manner  for  any  waveform  is  equal  to  the  rate  of  CO2  produc¬ 
tion  by  metabolism,  then  that  waveform  will  maintain  equilib¬ 
rium  of  arterial  CO2  partial  pressures.  Also  the  assumed 
boundary  condition  used  at  the  airway  opening  (i.e.,  PaoC02  = 
0  mmHg)  neglects  the  additional  dead  space  of  the  upper 
airways  or  an  endotracheal  tube.  Thus  we  expect  our  results  to 
slightly  underestimate  the  volume  amplitudes  required  for 


eucapnic  ventilation.  The  close  agreement  between  simulation 
results  and  the  experimental  data  of  Venegas  et  al.  (47)  may  be 
due  in  part  to  the  latter’s  use  of  specialized  oscillatory  appa¬ 
ratus  that  delivered  fresh  gas  directly  to  the  trachea,  bypassing 
the  dead  space  upper  airways  (44).  Boundary  conditions  at  the 
acinar  level  do  not  include  the  effects  of  acinar  interdepen¬ 
dence  (28),  collateral  ventilation  (4),  or  pleural  surface  inter¬ 
actions  at  the  chest  wall,  diaphragm,  or  between  lobes  (48). 

The  model  of  airway  segment  impedance  used  in  this  study 
(10)  was  intended  to  characterize  flow  phenomena  up  to  the 
lower  end  of  the  acoustic  range  (i.e.,  up  to  100  Hz),  assuming 
Womersley-type  oscillatory  flow  (42,  54)  and  isothermal  gas 
compression.  Admittedly,  the  lumped  longitudinal  and  shunt 
airway  segment  impedances  are  simplifications  of  actual  fluid 
mechanical  descriptions  of  oscillatory  flow.  The  isothermal  gas 
compression  compliance  used  in  our  model  is  also  a  limiting 
case  of  the  full  thermodynamic  expression  (6).  The  model  is 
also  limited  in  that  it  does  not  include  other  potentially  influ¬ 
ential  gas  transport  mechanisms,  such  as  asymmetric  velocity 
profiles  at  bifurcations,  differences  between  inspiratory  and 
expiratory  velocity  profiles,  or  cardiogenic  mixing  (9).  Further¬ 
more,  the  included  transport  mechanisms  are  assumed  to  su¬ 
perimpose  linearly,  except  for  direct  acinar  ventilation  (which 
is  inherently  nonlinear  due  to  complete  CO2  removal  at  the 
airway  opening).  Despite  limiting  the  scope  of  these  modeled 
mechanisms  for  the  sake  of  parsimony,  this  work  does  provide 
comparable  results  with  previous  experimental  work  using 
HFOV  in  dogs  (15,  27,  45). 

Aside  from  these  limitations,  there  are  several  important 
considerations  regarding  the  applicability  of  these  canine  sim¬ 
ulations  to  mechanical  ventilation  in  humans.  Differences  in 
anatomic  structure  of  the  canine  airway  tree  compared  with 
that  of  humans,  such  as  increased  branching  asymmetry  and 
dead  space  relative  to  body  weight,  will  also  yield  differences 
in  airway  resistance,  resonant  frequency,  and  the  transition 
frequency  from  conventional  ventilation  to  HFOV  (47).  More¬ 
over,  the  influence  of  frequency  on  the  pressure  cost  of  venti¬ 
lation,  as  well  as  on  the  distribution  of  gas  transport,  may  be 
altered  in  humans  compared  with  dogs  particularly  at  high 
frequencies  (46).  Clinically  relevant  predictions  for  optimal 
frequency  selection  during  HFOV  may  be  obtained  using  this 
gas  transport  model  in  airway  structures  based  on  human 
anatomy. 

Conclusion 

Oscillatory  frequency  has  a  substantial  impact  on  ventilation 
distribution.  The  total  tidal  volume  required  to  maintain  eu¬ 
capnia  decreases  with  frequency.  Thus  increasing  oscillatory 
frequency  may  minimize  the  risks  associated  with  excessive 
volume  distension.  However  increasing  frequency  may  not  be 
a  universally  appropriate  ventilator  management  strategy,  since 
frequency-dependent  transitions  from  elastically  dominated  to 
inertially  dominated  ventilation  distribution  results  in  overven¬ 
tilation  and  underventilation  simultaneously  in  different  lung 
regions.  Transitions  in  gas  transport  mechanisms  also  induce 
frequency-dependent  heterogeneity  in  ventilation  distribution, 
which  may  result  in  ventilation-to-perfusion  mismatching  and 
impairments  in  gas  exchange.  Thus  it  is  essential  to  consider 
these  multiple  interdependent  factors  when  selecting  frequency 
during  oscillation  of  the  injured  lung. 
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APPENDIX  A-1 

Derivation  of  mixing  transport  resistance 

Consider  an  example  bifurcation  given  by  three  cylindrical  airway 
segments  meeting  at  a  central  node.  The  CO2  concentration  at  the 
central  node  is  given  hy  C0CO2,  with  concentrations  at  the  distal  end 
of  each  of  the  three  adjoining  airway  segments  (i  =  1,  2,  3)  given  by 
C1CO2,  C2CO2,  and  C3CO2.  In  the  absence  of  diffusion  and  disper¬ 
sion,  the  time-averaged  value  of  C0CO2  is  given  by  a  weighted 
average  of  the  neighboring  concentrations  C,C02,  for  i  =  1,2,3.  Also 
assume  that  the  weighting  of  each  neighboring  CO2  concentration  is 
proportional  to  the  rate  at  which  the  mixing  volume  (Tmix,;)  passes 
completely  through  its  respective  airway  segment.  The  Kmix,!-  is  given 
hy  the  (strictly)  positive  difference  between  the  magnitude  of  total 
oscillating  volume  in  the  segment  IVoso,/l  and  the  segment  lumenal 
volume  itself  (Vseg,/): 

maX^  I  Vqsc,,' I  l^seg,/  ’  6)  (Ai) 


segment  (M““C02)  is  proportional  to  the  mixing  transport  resistance 
for  that  airway  The  axial  CO2  partial  pressure  difference  from 

the  distal  end  of  each  airway  segment  to  the  central  node  is  given  by 
PiC02  —  P0CO2.  Since  each  CO2  partial  pressure  differential  is 
defined  relative  to  the  central  node,  conservation  of  mass  of  the  CO2 
molar  flux  at  the  central  node  yields: 


2  Mf  "CO2 = i 

1=1  i=l 


P,C02  -  P0CO2 

nmix 

^T.i 


=  0 


(A4) 


Rearranging  Eq.  A4  to  solve  for  the  CO2  partial  pressure  at  the  central 
node  yields: 


/’0CO2  = 


(AS) 


where  max  ()  indicates  the  maximum  value  of  the  enclosed  arguments. 
The  weighted  averaged  concentration  at  the  central  node  is  then: 

i/V^mix.,C,.C02 


C0CO2  = 


i=l 


(A2) 


where /is  the  oscillation  frequency.  This  weighting  scheme  ensures 
that  i):  airway  segments  receiving  greater  proportions  of  incoming 
oscillatory  gas  flow  provide  greater  contributions  to  mixing  at  the 
bifurcation,  and  2)  airway  segments  receiving  oscillatory  gas  volumes 
less  than  their  respective  lumenal  volumes  do  not  contribute  to  mixing 
at  the  bifurcation.  Converting  CO2  concentrations  to  partial  pressures 
based  on  the  ideal  gas  law,  we  obtain: 


P0CO2 


,/V. 


mix,/ 


TZT 


P:C02 


jx  /Vmix.,- 


(A3) 


where  TZ  is  the  universal  gas  constant  and  T  is  the  gas  temperature.  We 
assume  that  this  mixing  transport  mechanism  can  be  characterized  by 
a  transport  resistance,  such  that  the  CO2  molar  flux  in  any  airway 


Combining  Eqs.  A3  and  A5,  it  can  be  demonstrated  that  Rjf  is 
inversely  proportional  to  the  weighting  factor: 

TZT 

=  (A6) 

J  ^mix,/ 


where  e  is  a  constant  of  proportionality  characterizing  the  complete¬ 
ness  with  which  Tmix,;  is  homogenized  between  neighboring  nodes. 
Very  large  s  yields  approaching  infinity,  resulting  in  no  mixing 
and  no  transport.  In  the  current  study  we  have  assumed  complete 
homogenization  of  Vmix.i  such  that  e  is  unity,  and  Eq.  A6  becomes: 


nmix 

^T,/ 


TZT 


fVm 


(AT) 


APPENDIX  A-2 

Mesh  analysis  of  transport  resistance  network 

The  transport  resistance  R™  of  an  airway  segment  relates  the  molar 
flux  M  and  CO2  partial  pressure  differential  ARCO2  in  the  axial 
direction.  For  a  network  of  branching  airways,  a  system  of  linear 
ordinary  differential  equations  can  be  easily  solved.  A  mesh  analysis 
of  the  system  ensures  conservation  of  mass,  such  that  the  sum  of 
molar  fluxes  through  the  airway  opening  and  every  terminal  airway 


Fig.  10.  Example  airway  network  and  corre¬ 
sponding  system  of  equations  for  solving 
molar  flux  of  CO2  out  of  each  acinus,  accord¬ 
ing  to  Eq.  A8.  This  network  comprises  5 
airways,  terminating  in  3  acini. 


tot 

T,5 
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segment  is  zero.  The  molar  fluxes  at  the  proximal  node  of  every  acinus 
are  chosen  as  state  variables,  and  the  system  of  equations  is  expressed 
in  matrix  notation  as: 

RjMC02  =  APCO2  (A8) 

where  each  row  of  the  system  corresponds  to  the  mesh  formed  hy  the 
airways  connecting  a  single  acinus  to  the  airway  opening.  Accord¬ 
ingly,  each  term  in  APCO2  corresponds  to  total  CO2  partial 
pressure  loss  between  an  nth  acinus  and  the  airway  opening  P„C02  — 
/'aoC02,  and  the  corresponding  nth  row  of  the  matrix  product 
ifxMCOa  corresponds  to  the  sum  of  partial  pressure  losses  due  to 
transport  resistance  in  each  airway  between  that  nth  acinus  and  the 
airway  opening. 

The  matrix  R-r  is  completely  dense,  because  the  transport  resistance 
of  the  trachea  f?x'trachea  included  in  every  mesh.  In  other  words, 
every  path  from  an  acinus  to  the  airway  opening  includes  the  same 
partial  pressure  loss  over  the  trachea,  which  is  equal  to 

V 

^T'tracheaS^;iC02  where  ti  indexes  all  individual  acini.  A  simple 

example  ol"  an  airway  network  and  corresponding  system  of  equations 
is  provided  in  Fig.  10.  Because  matrix  Rt  is  completely  dense, 
memory  cost  may  be  substantial  for  large  numbers  of  airways. 
However  Rt  is  symmetric  and  can  be  stored  and  solved  efficiently 
using  dedicated  algorithms  for  symmetric  matrices. 

An  alternative  to  mesh  analysis  is  nodal  analysis,  which  involves 
choosing  the  state  variables  as  the  partial  pressures  at  the  proximal 
node  of  every  airway  segment.  Nodal  analysis  forms  a  system  of 
equations  using  conservation  of  mass  applied  to  each  individual  node, 
which  produces  a  sparse  coefficient  matrix.  However,  small  numerical 
inaccuracies  at  each  node  can  result  in  substantial  discrepancies 
between  the  flux  at  the  airway  opening  and  the  total  flux  across  the 
acini.  Thus,  mesh  analysis  is  preferred. 

APPENDIX  A-3 

Expected  low-frequency  gas  exchange  behavior 

Expected  values  for  low-frequency  gas  exchange  (i.e.,  0.1  Hz  < 
/  <  1  Hz)  are  derived  from  standard  equations  of  gas  exchange 
physiology  during  conventional  mechanical  ventilation  and  spon¬ 
taneous  breathing.  Given  the  assumption  that  CO2  elimination  is 
dominated  by  direct  alveolar  ventilation  and  is  therefore  propor¬ 
tional  to  alveolar  ventilation,  we  have: 

VroTC02cc/(v^-- Vo)  (A9) 

Alveolar  partial  pressure  CO2  is  related  to  the  eucapnic  alveolar 

ventilation  (Va)  and  CO2  elimination  rate  (V^“C02)  and  is  assumed  to 
be  equal  to  mixed  arterial  CO2  partial  pressure  (PaC02): 

1^“C02  V:“C02  7’aC02 

Va  863 

where  863  is  a  conversion  factor  between  partial  pressure  in  mmHg 
and  gas  volume  fraction,  accounting  for  the  difference  between 

C=C02  in  1/s  at  standard  temperature  and  pressure  dry  air,  and  Va  in 
1/s  at  body  temperature  and  pressure  air  saturated  with  water  vapor. 
Equation  AlO  can  be  rearranged  in  the  form  of  the  power-law 
regression  in  Eq.  15,  such  that: 

/  863  .  \  , 

+  (  ^:^^oTC02 jr '  (All) 

Thus,  the  expected  values  for  the  power-law  regression  parameters  for 
/  863  .  \ 

Eq.  15  are  Po  =  Vd,  Pi  =  I  ^^C'C02j,  and  P2  =  -1. 


ACKNOWLEDGMENTS 

Portions  of  this  work  were  excerpted  from  Master  of  Science  thesis 
“Simulation  of  Ventilation  Distribution  and  Gas  Transport  During  Oscillatoiy 
Ventilation”  by  Jacob  Herrmann,  Boston  University,  2015.  The  authors  thank 
Drs.  Bela  Suki,  Brett  Simon,  J.  Jane  Pillow,  and  Joseph  Tien  for  helpful 
suggestions  during  the  course  of  this  work. 

GRANTS 

This  work  was  supported  in  part  by  the  Medical  Technologies  Centre  of 
Research  Excellence  at  the  University  of  Auckland  (M.  H.  Tawhai),  National 
Institutes  of  Health  Grant  UMl  HL-108724  (D.  W.  Kaczka),  US  Department 
of  Defense  Grant  PR151761  (D.  W.  Kaczka),  and  the  University  of  Iowa 
Department  of  Anesthesia  (D.  W.  Kaczka,  J.  Herrmann). 

DISCLOSURES 

D.  W.  Kaczka  and  J.  Heixmann  are  coinventors  on  a  pending  patent 
involving  multifrequency  oscillatory  ventilation  (MFOV).  In  addition,  they  are 
cofounders  and  shareholders  of  Oscillavent,  Inc. 

AUTHOR  CONTRIBUTIONS 

J.H.,  M.H.T.,  and  D.W.K.  conception  and  design  of  research;  J.H.  and 
M.H.T.  performed  experiments;  J.H.,  M.H.T.,  and  D.W.K.  analyzed  data;  J.H., 
M.H.T.,  and  D.W.K.  interpreted  results  of  experiments;  J.H.  and  D.W.K. 
prepared  figures;  J.H.  and  D.W.K.  drafted  manuscript;  J.H.,  M.H.T.,  and 
D.W.K.  edited  and  revised  manuscript;  J.H.,  M.H.T.,  and  D.W.K.  approved 
final  version  of  manuscript. 

REFERENCES 

1.  The  Acute  Respiratory  Distress  Syndrome  Network.  Ventilation  with 
lower  tidal  volumes  as  compared  with  traditional  tidal  volumes  for  acute 
lung  injury  and  the  acute  respiratory  distress  syndrome.  N  Engl  J  Med  342: 
1301-1308,  2000.  doi:10.1056/NEJM200005043421801. 

2.  Amini  R,  Herrmann  J,  Kaczka  DW.  Intratidal  overdistention  and  dere¬ 
cruitment  in  the  injured  lung:  A  simulation  study  [epub  ahead  of  print]. 
IEEE  Trans  Biomed  Eng  in  press,  2016.  doi:  10.1 109/TBME. 
2016.2572678. 

3.  Amini  R,  Kaczka  DW.  Impact  of  ventilation  frequency  and  parenchymal 
stiffness  on  flow  and  pressure  distribution  in  a  canine  lung  model.  Ann 
Biomed  Eng  41:  2699-2711,  2013.  doi:10.1007/sl0439-013-0866-7. 

4.  Armengol  J,  Jones  RL,  King  EG.  Collateral  ventilation  during  high- 
frequency  oscillation  in  dogs  J  Appl  Physiol  (1985)  58:  173-179,  1985. 

5.  Bates  JHT,  Irvin  CG.  Time  dependence  of  recruitment  and  derecruitment 
in  the  lung:  a  theoretical  model.  J  Appl  Physiol  ( 1985}  93:  705-713,  2002. 
doi:  10.11 52/japplphysiol.0 1 274.200 1 . 

6.  Benade  AH.  On  the  propagation  of  sound  waves  in  a  cylindrical  conduit. 
JAcoust  SocAm  44:  616-623,  1968.  doi:10.1121/l. 1911130. 

7.  Berdine  GG,  Lehr  JL,  McKinley  DS,  Drazen  JM.  Nonuniformity  of 
canine  lung  washout  by  high-frequency  ventilation.  J  Appl  Physiol  (1985) 
61:  1388-1394,  1986. 

8.  Carvalho  AR,  Spieth  PM,  Pelosi  P,  Vidal  Melo  MF,  Koch  T,  Jandre 
FC,  Giannella-Neto  A,  de  Abreu  MG.  Ability  of  dynamic  airway 
pressure  curve  profile  and  elastance  for  positive  end-expiratory  pressure 
titration.  Intensive  Care  Med  34:  2291-2299,  2008.  doi:10.1007/s00134- 
008-1301-7. 

9.  Chang  HK.  Mechanisms  of  gas  transport  during  ventilation  by  high- 
frequency  oscillation.  J  Appl  Physiol  Respir  Environ  Exerc  Physiol  56: 
553-563,  1984. 

10.  Colletti  AA,  Amini  R,  Kaczka  DW.  Simulating  ventilation  distribution  in 
heterogenous  lung  injury  using  a  binary  tree  data  structure.  Comput  Biol 
Med  41:  936-945,  2011.  doi:10.1016/j.compbiomed.2011.08.004. 

11.  Ebert  DS,  Musgrave  FK,  Peachey  D,  Perlin  K,  Worley  S,  Mark  WR, 
Hart  JC.  Texture  and  Modeling:  A  Procedural  Approach  (3rd  ed.)-  San 
Francisco,  CA:  Morgan  Kaufmann  Publishers,  2002. 

12.  Ferguson  ND,  Cook  DJ,  Guyatt  GH,  Mehta  S,  Hand  L,  Austin  P,  Zhou 
Q,  Matte  A,  Walter  SD,  Lamontagne  F,  Granton  JT,  Arabi  YM, 
Arroliga  AC,  Stewart  TE,  Slutsky  AS,  Meade  MO;  OSCILLATE 
Trial  Investigators;  Canadian  Critical  Care  Trials  Group.  High- 
frequency  oscillation  in  early  acute  respiratory  distress  syndrome.  N  Engl 
J  Med  368:  795-805,  2013.  doi:10.1056/NEJMoal215554. 

13.  Fortune  JB,  Wagner  PD.  Effects  of  common  dead  space  on  inert  gas 
exchange  in  mathematical  models  of  the  lung.  J  Appl  Physiol  Respir 
Environ  Exerc  Physiol  47:  896-906,  1979. 


J  Appl  Physiol  •  doi:  10. 1 152/japplphysiol. 00097. 2016  •  www.jappl.org 


Downloaded  from  http://jap.physiology.org/  by  10.220.32.246  on  December  23,  2016 


1318 


Gas  Transport  during  Oscillatory  Ventilation  •  Hemnann  J  et  al. 


14.  Fredberg  JJ.  Augmented  diffusion  in  the  airways  can  support  pulmonaiy 
gas  exchange.  J  Appl  Physiol  Respir  Environ  Exerc  Physiol  49:  232-238, 
1980. 

15.  Fredberg  JJ,  Keefe  DH,  Glass  GM,  Castile  RG,  Frantz  ID  III. 
Alveolar  pressure  nonhomogeneity  during  small-amplitude  high-fre¬ 
quency  oscillation.  J  Appl  Physiol  Respir  Environ  Exerc  Physiol  57: 
788-800,  1984. 

16.  Greenblatt  EE,  Butler  JP,  Venegas  JG,  Winkler  T.  Pendelluft  in  the 
bronchial  tree.  J  Appl  Physiol  (1985)  117:  979-988,  2014.  doi:10.1152/ 
japplphysiol.00466.2014. 

17.  Hantos  Z,  Daroczy  B,  Suki  B,  Nagy  S,  Fredberg  JJ.  Input  impedance 
and  peripheral  inhomogeneity  of  dog  lungs.  J  Appl  Physiol  (1985)  72: 
168-178,  1992. 

18.  Haskins  S,  Pascoe  PJ,  Ilkiw  JE,  Fudge  J,  Hopper  K,  Aldrich  J. 

Reference  cardiopulmonary  values  in  normal  dogs.  Comp  Med  55:  156- 
161,  2005. 

19.  Jaeger  MJ.  High-frequency  ventilation  in  dogs  with  three  gases  of 
different  densities.  J  Appl  Physiol  (1985)  70:  2188-2192,  1991. 

20.  Jaeger  MJ,  Kurzweg  UH,  Banner  MJ.  Transport  of  gases  in  high- 
frequency  ventilation.  Crit  Care  Med  12:  708—710,  1984.  doi:10.1097/ 
00003246-198409000-00003. 

21 .  Johansen  T,  Winkler  T,  Kelly  VJ,  Osorio-Valencia  JS,  Greenblatt  EE, 
Harris  RS,  Venegas  JG.  A  method  for  mapping  regional  oxygen  and 
C02  transfer  in  the  lung.  Respir  Physiol  Neurobiol  222:  29-47,  2016. 
doi:10.1016/j.resp.2015. 10.017. 

22.  Kaczka  DW,  Brown  RH,  Mitzner  W.  Assessment  of  heterogeneous 
airway  constriction  in  dogs:  a  structure-function  analysis.  J  Appl  Physiol 
(1985)  106:  520-530,  2009.  doi:10.1152/japplphysiol.90576.2008. 

23.  Kaczka  DW,  Herrmann  J,  Zonneveld  CE,  Tingay  DG,  Lavizzari  A, 
Noble  PB,  Pillow  JJ.  Multifrequency  oscillatory  ventilation  in  the 
premature  lung:  effects  on  gas  exchange,  mechanics,  and  ventilation 
distribution.  Anesthesiology  123:  1394-1403,  2015.  doi:10.1097/ 
ALN.0000000000000898. 

24.  Kaczka  DW,  Massa  CB,  Simon  BA.  Reliability  of  estimating  stochastic 
lung  tissue  heterogeneity  from  pulmonary  impedance  spectra:  a  forward- 
inverse  modeling  study.  Ann  Biomed  Eng  35:  1722-1738,  2007.  doi: 
1 0. 1 007/s  1 0439-007-9339- 1 . 

25.  Khoo  MC,  Slutsky  AS,  Drazen  JM,  Solway  J,  Gavriely  N,  Kamm  RD. 
Gas  mixing  during  high-frequency  ventilation:  an  improved  model.  J  Appl 
Physiol  Respir  Environ  Exerc  Physiol  57:  493-506,  1984. 

26.  Lai-Fook  SJ,  Rodarte  JR.  Pleural  pressure  distribution  and  its  relation¬ 
ship  to  lung  volume  and  interstitial  pressure.  J  Appl  Physiol  (1985)  70: 
967-978,  1991. 

27.  Lehr  JL,  Butler  JP,  Westerman  PA,  Zatz  SL,  Drazen  JM.  Photo¬ 
graphic  measurement  of  pleural  surface  motion  during  lung  oscillation.  J 
Appl  Physiol  (1985)  59:  623-633,  1985. 

28.  Mead  J,  Takishima  T,  Leith  D.  Stress  distribution  in  lungs:  a  model  of 
pulmonary  elasticity.  J  Appl  Physiol  28:  596-608,  1970. 

29.  Otis  AB,  McKerrow  CB,  Bartlett  RA,  Mead  J,  Mcllroy  MB,  Selver- 
Stone  NJ,  Radford  EP  Jr.  Mechanical  factors  in  distribution  of  pulmo¬ 
nary  ventilation  [Online].  [8  Sep.  2014].  J  Appl  Physiol  8:  427—443, 
1956.http://www.ncbi.nlm.nih.gov/pubmed/l  3286206 

30.  Piiper  J,  Scheid  P.  Diffusion  and  convection  in  intrapulmonaiy  gas 
mixing,  in  Comprehensive  Physiology.  Bethesda,  MD:  John  Wiley  & 
Sons,  p.  51-69.  doi:10.1002/cphy.cp030404. 

31.  Pillow  JJ.  High-frequency  oscillatory  ventilation:  mechanisms  of  gas 
exchange  and  lung  mechanics.  Crit  Care  Med  33,  Suppl:  S135-S141, 
2005.  doi:10.1097/01.CCM.0000155789.52984.B7. 

32.  Poulsen  P,  Karbing  DS,  Rees  SE,  Andreassen  S.  Tidal  breathing  model 
describing  end-tidal,  alveolar,  arterial  and  mixed  venous  CO2  and  O2. 
Comput  Methods  Programs  Biomed  101:  166-172,  2011.  doi:10.1016/ 
j.cmpb.2010.03.020. 

33.  Scherer  PW,  Haselton  FR,  Seybert  JR.  Gas  transport  in  branched 
airways  during  high-frequency  ventilation.  Ann  Biomed  Eng  12:  385^05, 
1984.  doi:10.1007/BF02407782. 

34.  Scherer  PW,  Shendalman  LH,  Greene  NM,  Bouhuys  A.  Measurement 
of  axial  diffusivities  in  a  model  of  the  bronchial  airways.  J  Appl  Physiol 
38:  719-723,  1975. 


35.  Slutsky  AS,  Drazen  FM,  Ingram  RH  Jr,  Kamm  RD,  Shapiro  AH, 
Fredherg  JJ,  Loring  SH,  Lehr  J.  Effective  pulmonary  ventilation  with 
small-volume  oscillations  at  high  frequency.  Science  209:  609-671,  1980. 
doi:  10.11 26/science. 677 1872. 

36.  Slutsky  AS,  Kamm  RD,  Rossing  TH,  Loring  SH,  Lehr  J,  Shapiro  AH, 
Ingram  RH  Jr,  Drazen  JM.  Effects  of  frequency,  tidal  volume,  and  lung 
volume  on  C02  elimination  in  dogs  by  high  frequency  (2-30  Hz),  low  tidal 
volume  ventilation.  J  Clin  Invest  68:  1475-1484,  1981.  doi:10.1172/ 
JCl  110400. 

37.  Slutsky  AS,  Ranieri  VM.  Ventilator-induced  lung  injury.  N  Engl  J  Med 
369:  2126-2136,  2013.  doi:10.1056/NEJMral208707. 

38.  Tawhai  MH,  Hunter  P,  Tschirren  J,  Reinhardt  J,  McLennan  G, 
Hoffman  EA.  CT-based  geometry  analysis  and  finite  element  models  of 
the  human  and  ovine  bronchial  tree.  J  Appl  Physiol  ( 1985)  97:  2310-2321, 
2004.  doi:  10. 1 1 52/japplphysiol. 00520. 2004. 

39.  Howatson  Tawhai  M,  Pullan  AJ,  Hunter  PJ.  Generation  of  an  anatom¬ 
ically  based  three-dimensional  model  of  the  conducting  airways.  Ann 
Biomed  Eng  28:  793-802,  2000.  doi:10.1114/l. 1289457. 

40.  Taylor  G.  The  dispersion  of  matter  in  turbulent  flow  through  a  pipe.  Proc 
R  Soc  A  Math  Phys.  Eng  Sci  223:  446^68,  1954. 

41.  Thorpe  CW,  Bates  JHT.  Effect  of  stochastic  heterogeneity  on  lung 
impedance  during  acute  bronchoconstriction:  a  model  analysis.  J  Appl 
Physiol  (1985)  82:  1616-1625,  1997. 

42.  Thurston  GB.  Periodic  fluid  flow  through  circular  tubes.  J  Acoust  Soc  Am 
24:  653-656,  1952.  doi:10.1 121/1.1906949. 

43.  Tsuzaki  K,  Hales  CA,  Strieder  DJ,  Venegas  JG.  Regional  lung  me¬ 
chanics  and  gas  transport  in  lungs  with  inhomogeneous  compliance.  J 
Appl  Physiol  (1985)  75:  206-216,  1993. 

44.  Venegas  JG.  Equivalent  circuit  analysis  of  high-frequency  ventilators 
including  a  new  high-impedance  flow-interrupting  ventilator.  IEEE  Trans 
Biomed  Eng  33:  420-427,  1986.  doi:10.1109/TBME.1986.325798. 

45.  Venegas  JG,  Custer  J,  Kamm  RD,  Hales  CA.  Relationship  for  gas 
transport  during  high-frequency  ventilation  in  dogs.  J  Appl  Physiol  (1985) 
59:  1539-1547,  1985. 

46.  Venegas  JG,  Fredberg  JJ.  Understanding  the  pressure  cost  of  ventilation: 
why  does  high-frequency  ventilation  work?  Crit  Care  Med  22,  Supph 
S49-S57,  1994.  doi:10.1097/00003246-199422091-00004. 

47.  Venegas  JG,  Hales  CA,  Strieder  DJ.  A  general  dimensionless  equation 
of  gas  transport  by  high-frequency  ventilation.  J  Appl  Physiol  (1985)  60: 
1025-1030,  1986. 

48.  Venegas  JG,  Tsuzaki  K,  Fox  BJ,  Simon  BA,  Hales  CA.  Regional 
coupling  between  chest  wall  and  lung  expansion  during  HFV :  a  positron 
imaging  .study.  J  Appl  Physiol  (1985)  1A\  2242-2252,  1993. 

49.  Venegas  JG,  Yamada  Y,  Custer  J,  Hales  CA.  Effects  of  respiratory 
variables  on  regional  gas  transport  during  high-frequency  ventilation.  J 
Appl  Physiol  (1985)  64:  2108-2118,  1988. 

50.  Watson  EJ.  Diffusion  in  oscillatory  pipe  flow.  J  Fluid  Mech  133: 
233-244,  1983.  doi:10.1017/S0022112083001883. 

51.  Watson  JW,  Jackson  AC.  C02  elimination  as  a  function  of  frequency 
and  tidal  volume  in  rabbits  during  HFO.  J  Appl  Physiol  Respir  Environ 
Exerc  Physiol  57:  354-359,  1984. 

52.  Watson  JW,  Jackson  AC.  Frequency  dependence  of  C02  elimination 
and  respiratory  resistance  in  monkeys.  J  Appl  Physiol  (1985)  58:  653-657, 
1985. 

53.  West  JB.  Ventilation-perfusion  inequality  and  overall  gas  exchange  in 
computer  models  of  the  lung.  Respir  Physiol  7:  88-110,  1969.  doi: 
10.1016/0034-5687(69)90071-1. 

54.  Womersley  JR.  Method  for  the  calculation  of  velocity,  rate  of  flow  and 
viscous  drag  in  arteries  when  the  pressure  gradient  is  known.  J  Physiol 
127:  553-563,  1955.  doi:10.1113/jphysiol.l955.sp005276. 

55.  Yamada  Y,  Burnham  C,  Hales  CA,  Venegas  JG.  Regional  mapping  of 
gas  transport  during  high-frequency  and  conventional  ventilation.  J  Appl 
Physiol  (1985)  66:  1209-1218,  1989. 

56.  Young  D,  Lamb  SE,  Shah  S,  MacKenzie  I,  Tunnicliffe  W,  Lall  R, 
Rowan  K,  Cuthbertson  BH;  OSCAR  Study  Group.  High-frequency 
oscillation  for  acute  respiratory  distress  syndrome.  N  Engl  J  Med  368: 
806-813,  2013.  doi:10.1056/NEJMoal215716. 


J  Appl  Physiol  •  doi:  10. 1 1 52/japplphysiol. 00097. 2016  •  www.jappl.org 


Downloaded  from  http://jap.physiology.org/  by  10.220.32.246  on  December  23,  2016 


4/3/2017 


cOASIS,  The  Online  Abstract  Submission  System 


2017  Meeting  of  the  American  Society  of  Anesthesiologists 

ANESTHESSOLOGY 

BOSTON 


lO 


OCTOBER  2 1-25 


Print  this  Page  for  Your  Records 


Ciose  Window 


Control/Tracking  Number  :  17-SA-4940-ASAHQ 
Activity  iScientific  Abstract 
Current  Date/Time  :  4/3/2017  8:09:43  PM 

TiTLE: 

Parenchymai  Strain  vs.  Strain  Rate  During  Muiti-Frequency  Osciiiatory  Ventiiation 

AUTHOR(S): 

David  W.  Kaczka,  M.D.,Ph.D.\  Jacob  Herrmann,  M.S.^,  Merryn  H.  Tawhai,  Ph.D.^. 

^Anesthesia,  Biomedicai  Engineering,  and  Radioiogy,  University  of  Iowa,  Iowa  City,  lA,  USA,  ^Anesthesia  and  Biomedicai  Engineering,  University 
of  Iowa,  Iowa  City,  lA,  USA,  ^Auckiand  Bioengineering  Institute,  University  of  Auckiand,  Auckiand,  New  Zeaiand. 


AFFIRMATIONS: 

Affirmations  Cont.  (Complete): 

I  Agree:  True 

*:  I  agree  to  the  above  statements 
*:  No  human  subjects  were  invoived  in  the  research 
*:  No  animai  subjects  were  invoived  in  the  research 
*:  No  human  or  animai  subjects  were  invoived  in  the  research 


SESSION  CATEGORY: 

017.2  RESPIRATION  -  Pulmonary  Mechanics  and  Ventilation  (Including  Artificial) 


QUESTIONNAIRE: 

Questionnaire  (Complete): 

Please  select:  No,  do  not  consider  my  abstract  for  the  ANESTHESIOLOGY  Journal  Symposium 

Is  this  abstract  submitted  by  a  resident?:  No 
Funding  Information:  Other  Federal  Funding 

If  you  chose  "Other"  above,  please  explain  (250  characters):  U  S.  Department  of  Defense 

Yes  or  No?:  No 

Please  select:  Yes,  I  would  like  to  be  considered  to  present  an  Oral/Poster  Discussion  Session. 

Please  Select:  If  my  abstract  is  selected  for  a  poster  presentation  at  ANESTHESIOLOGY*  2017,  I  would  have  no  preference  between  an 

electronic  or  a  printed  poster  presentation 


ABSTRACT: 

Introduction:  Ventilator-induced  lung  injury  (VILI)  in  mechanically  ventilated  patients  may  result  from  extreme  stresses,  strains,  and  strain  rates 
associated  with  parenchymal  overdistension,  as  well  as  repetitive  alveolar  collapse  and  reopening.  High  frequency  oscillatory  ventilation  (HFOV) 
is  an  alternative  lung-protective  ventilation  strategy,  which  attempts  to  minimize  VILI  using  small  volume  amplitudes  delivered  at  high 
frequencies  (2  to  20  Hz).  However  oscillatory  ventilation  at  a  single  discrete  frequency  is  distributed  non-uniformly  and  in  a  frequency- 
dependent  manner,  especially  in  the  mechanically  heterogeneous  lung.  In  a  previous  study,  we  demonstrated  that  lung  function  and  gas 
exchange  can  be  significantly  improved  if  volume  oscillations  are  applied  at  multiple  simultaneous  frequencies,  rather  than  at  a  single  discrete 
frequency  (Kaczka  et  al.  Anesthesiology  123:1394-1403,  2015).  We  refer  to  this  unique  ventilatory  modality  as  ‘Multi-Frequency  Oscillatory 
Ventilation'  (MFOV),  and  propose  that  improved  physiologic  outcomes  arise  from  the  more  even  distribution  of  ventilation,  in  accordance  with 
local  lung  mechanical  properties.  In  the  present  study,  we  hypothesized  that  the  spectral  content  of  MFOV  waveforms  can  be  further  tuned 
according  to  the  relative  contributions  of  acinar  strain  vs.  strain  rate  to  VILI.  We  developed  a  computational  model  of  the  porcine  lung  to 
investigate  the  potential  for  minimizing  distributed  acinar  strains  and  strain  rates  during  MFOV.  Methods:  The  computational  model  was 
constructed  using  a  central  airway  tree  segmented  from  a  computed  tomographic  image  of  a  supine  pig,  with  algorithmically  generated 
peripheral  airways.  The  model  consisted  of  60,494  airway  segments  in  total,  with  30,243  viscoelastic  acini.  Eucapnic  ventilation  was  simulated 
using  a  Monte  Carlo  technique  to  select  from  over  one  thousand  MFOV  waveforms  consisting  of  four  frequencies  (5,  10,  15,  and  20  Hz),  with 
randomized  volume  amplitudes  and  phases.  Corresponding  distributions  of  acinar  strains  and  strain  rates  throughout  the  model  were 
determined  for  each  MFOV  waveform,  and  used  to  compute  a  cost  function  defined  by  weighting  terms  to  emphasize  the  percent  contributions 
of  parenchymal  strain  vs.  strain  rate  to  VILI.  Results:  Our  simulations  demonstrated  that  MFOV  waveforms  were  superior  to  traditional  single¬ 
frequency  HFOV  for  minimizing  strain  and  strain  rate  in  both  healthy  and  injured  lungs.  Optimized  volume  amplitudes  at  each  frequency  in  the 
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MFOV  waveform  were  modulated  according  to  the  relative  cost  function  weighting  of  strain  and  strain  rate,  as  well  as  the  degree  of  lung 
mechanical  heterogeneity  (Figure  1).  Increasing  contribution  of  strain  rate  to  VILI  (as  opposed  to  strain)  resulted  in  decreasing  volume 
amplitudes  at  higher  frequencies  for  the  optimized  MFOV  waveforms.  Conclusion:  These  simulations  indicate  that  frequency  content  in  MFOV 
waveforms  may  be  tuned  to  minimize  the  relative  contributions  of  parenchymal  strain  and  strain  rate.  Our  unique  approach  may  allow  for  the 
selection  of  patient-specific  optimized  MFOV  waveforms,  especially  when  combined  with  experimental  evidence  to  justify  physiologically- 


relevant  emphasis  on  strain  vs.  strain  rate  to  minimize  risk  for  VILI. 


Cost  function  weig^tlngon  strain  (%) 


Figure  1:  Computational  model  of  porcine  lung  (left  panel).  Optimal  volume  amplitudes  at  each 
frequency  are  shown  for  varying  cost  functon  emphasis  on  strain  vs.  strain  rate,  in  simulated 
healthy  conditions  (center  panel)  and  heterogeneous  lung  injury  (right  panel). 
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SUMMARY: 

Multi-Frequency  Oscillatory  Ventilation  (MFOV)  has  been  demonstrated  to  improve  physiologic  outcomes  in  heterogeneous  lungs,  due  to  more 
even  distribution  of  ventilation.  In  this  computational  modeling  study  of  the  porcine  lung,  we  demonstrate  that  the  spectral  content  of  MFOV 
waveforms  can  be  further  optimized,  according  to  the  relative  contributions  of  parenchymal  strain  vs.  strain  rate  to  ventilator-induced  lung 

injury. 
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